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A COMPUTER PROGRAM FOR CALCULATING LAMINAR AND TURBULENT 
BOUNDARY LAYERS FOR TWO-DIMENSIONAL TIME-DEPENDENT FLOWS 

by 

Tuncer Cebeci* and Lawrence W. Carr 


SUMMARY 

There are many unsteady boundary-layer problems in which it is necessary 
to account for fluctuations in the external flow. These fluctuations may 
change in direction and magnitude and, in most turbulent flow cases, may be 
regarded as low frequency fluctuations superimposed on the turbulence energy 
spectrum. The computer program described here provides solutions of two- 
dimensional equations appropriate to laminar and turbulent boundary layers for 
boundary conditions with an external flow which fluctuates in magnitude; it 
can readily be extended to the more general case. It is based on the numerical 
solution of the governing boundary-layer equations by an efficient two-point 
finite-difference method developed by Keller and Cebeci [1]. An eddy- 
viscosity formulation is used to model the Reynolds shear-stress term. 

Here, after briefly describing the main features of the method, we provide 
instructions for the computer program with a listing, and present sample 
calculations to demonstrate its usage and capabilities for laminar and 
turbulent unsteady boundary layers with an external flow which fluctuates in 
magnitude. For further details of the metnod, the reader is referred to 
reference 2. 


California State University, 1250 N. Bellflower Blvd., 

Long Beach, Calif. 90840. 
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I. INTRODUCTION 


The calculation of the properties of unsteady boundary layers is important 
to a wide range of applications including the blades of compressors, turbines 
and helicopters. In each of these cases, the wakes of previous blades generate 
a freestream velocity for a subsequent blade /ith regular, time-dependent fluctu- 
ations in amplitude and perhaps frequency. In the helicopter arrangement, 
variable pitch also leads to a phase relationship between the blade surface and 
the freestream flow. Discrete-frequency fluctuations of this type influence the 
characteristics of boundary-layer flows at all Reynolds numbers. The present 
report is concerned with a method, previously used by Cebeci [2], to allow the 
efficient and accurate calculation of two-dimensional, boundary-layer flows in 
the presence of a freestream which varies in amplitude at uniform frequencies. 

The calculation method is, however, general and is readily capable of extension 
to include freestream velocity distributions which vary in amplitude and fre- 
quency, albeit over a limited range of frequencies; it can also be extended to 
represent the time-dependent flow over arbitrary three-dimensional bodies. 

The solution of the two-dimensional, boundary-layer equations with constant 
viscosity, an unsteady convective term and a prescribed and unsteady longitudinal 
pressure gradient is not novel and corresponding calculations are included here 
to allow the capabilities of the present numerical procedure to be compared with 
alternatives. McCroskey [3], in a recent and comprehensive review, has suggested 
that the earlier reviews of Riley [4] and Patel [5] provide an essentially com- 
plete theoretical understanding of unsteady laminar flow in the absence of strong 
pressure gradients. The present method can be used to generate solutions for 
strong pressure-gradient situations, if required and, with the aid of the Mechul- 
function approach of reference 6, to calculate through small regions of separated 
flow. In flows where the direction of velocity reverses and solutions can be 
obtained without inverse methods, the zig-zag modification to the Box scheme 
used by Cebeci [7] to calculate the properties of the laminar flow over a cylinder 
started impulsively from rest, is likely to be necessary. 

An important emphasis of the present work is on the application of an effici- 
ent, accurate and general numerical procedure to the calculation of unsteady 
boundary- layer flows. Accuracy, with efficiency, is necessary if practical 





calculations are to be performed through regions of flow separation or over three- 
dimensional geometries and the present two-dimensional calculations provide a 
basis for the mor.- general purposes. 

In formulating equations for turbulent flow which are similar in form to 
those for laminar flow, it is presumed that the eddy viscosity represents all 
turbulent viscous effects. The unsteady effects of the freestream are super- 
imposed on the turbulence spectrum through the time-dependent convective term. 
Alternative, but invariably more complicated, approaches are possible. For 
example, turbulence energy equations can be solved with frequency or wave number 
as dependent variables or the subgrid-scale modelling approach can be used with 
solutions of time-dependent equations only at the lower frequencies. The present 
eddy-viscosity formulation has the advantage of simplicity and has a wide range 
of applicability. It is likely, however, that hiqher-order turbulence models 
will be required to represent flows with strong favorable and adverse pressure 
gradients even though the present method can be adequate for some separating 
boundary layers as shown here. 

The report has been prepared with three main sections describing, respec- 
tively, the calculation method, the computer program and sample calculations. 

The description of the calculation method is brief since the equations, the 
turbulence model, coordinate transformations and solution procedure have been 
discussed in detail, for example in reference 2. The computer program is 
described in sufficient detail to allow it to be used and the samp.e calcula- 
tions are also described in a manner which is intended to assist a user. The 
results presented, as a consequence of the sample calculations, demonstrate 
that the program can conveniently lead to precise results for steady and unsteady 
laminar and turbulent flows and suggest that future extensions to include free- 
stream velocities which vary in amplitude and frequency, to separated flow and 
to three-dimensional flow can readily be achieved. There is, however, an inade- 
quate range of experimental data against which to evaluate the calculations. 

The report ends with a brief discussion and summary conclusions. 
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II. DESCRIPTION OF THE METHOD 


2.1 Boundary-Layer Equations in Physical Variables 

The governing boundary-layer equations for incompressible unsteady 
laminar and turbulent boundary layers are: 
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and are considered here with the following boundary conditions for the wall 
and for the outer edge of the boundary layer: 
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To complete the formulation of the problem, initial conditions must be 
specified in the (t,y) plane for x = x Q and in the (x,y) plane for t = t Q . 
Here we consider a flow in which at time t = 0, the flow field is given by 
steady-state conditions; and, for t > 0, the external velocity u g (x,t) 
fluctuates from the steady-state velocity u Q (x) according to the expression 

u e i x , t ) « u 0 (x)(l + B coswt) (4) 

At time t = 0, the flow can be either laminar or turbulent. In the latter 
case, the surface distance x must be greater than zero. The initial con- 
ditions in the (t,y) plane at x = 0 are obtained from the similarity 
solutions of (1) to (3) for laminar flows as described in detail in 
reference 2 and later in this report. 


2.2 Turbulence Model 

The present calculations use the eddy-viscosity concept to model the Reynolds 
shear stress term of equation (2). According to the formulation described in 
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reference [2], the eddy viscosity is defined by 

-t — r 3u 
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with two separate formulas used to represent the inner and outer regions of the 
boundary layer. In the so-called inner region of the boundary layer, e m is 
defined by the formula 

(e m ) i - {0.4y [1 -exp(-y/A)]} 2 N (6) 


where 
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and z-j = Rq/425 - 1 

The inner and outer regions are established by the continuity of the eddy- 
viscosity formulas. 

2.3 Governing Equations in Transformed Variables 


Before the governing system of equations described above is solved by 
the numerical method to be described in the following section, they are first 
expressed in transformed variables which are defined by 


x = x. 


t = t, 
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In addition, a dimensionless stream function f(x,t,n) defined by 

'P 88 (vxu Q )^ 2 f(x,t,n) 

is used to satisfy the continuity equation and to express the momentum equa- 
tion (2) as a third-order equation. Here u Q Is some reference velocity 
and ij/ is the stream function defined by 
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With equations (10) and (11) and with the definition of eddy viscosity, the 
continuity and momentum equations, (1) to (2), and their boundary conditions 
(3) may be written as 
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Here primes denote differentiation with respect to n and 
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The initial conditions at x = 0 in the (n, t) plane for a laminar flow 
are obtained by writing (12) as 
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Note that, although for a flat-plate flow the right-hand side of (15) is 
equal to zero, for stagnation-point flow (u Q = Ax) it is not. For this 
reason, the slope A(^du g /dx) must be specified for flows which start with 
a stagnation point. 


The initial conditions at t = 0 in the (n.x)-plane for a laminar or 
turbulent flow are obtained by writing equation (12) as 


(bf")' + ff" -P(f ') 2 + P 3 = x(f f£--f" f£y 


(16) 





Equation (16) expresses conservation of momentum in transformed variables for 
steady-state conditions. Here P 3 is given by 
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In terms of transformed variables, the eddy viscosity formulas become 
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where the subscript « refers to the boundary-layer edge and the parameters 
R x and y/A are given by 

u.x Rl / 4 (f ;) 1/2 . . 1/2 

R x = -v- & = J Sr — "[i - n.8(p+ + p*)3 1/2 09) 


2.4 Solution Procedure 

We use the numerical method of reference 1 to solve the system of equa- 
tions described in the previous sections. This is an efficient two-point 
finite-difference method developed by Keller and Cebeci tH and extensively used 
by Cebeci for two-dimensional and three-dimensional flows (see for example, 
references 8 and 9). A detailed description is presented in references 8 and 
10 and is not repeated here. 

One of the advantages of the numerical method 4 s that nonuniform net 
spacings can be used in the t, x directions as well as across the boundary 
layer. In the latter case, the nonuniform grid is a geometric progression 
with the property that the ratio of lengths of any two adjacent Intervals is 
a constant; that is, Anj = KArij_-j - The distance to the j-th line is given 
by the following formula: 

An j - Arv,(K J - 1)/(K - 1) K > 1 (20) 

There are two parameters in equation (20): Ar^ , the lenqth of the first step, 

and K, the ratio of two successive steps. The total number of points J can 
be calculated from the following formula: 
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In the computer program which embodies the present solution method, Ar^ and 
K are chosen with typical values, for moderate Reynolds numbers, of 0.01 
and 1.3, respectively. In general, approximately 50 grid nodes across the 
boundary layer are sufficient to represent laminar and turbulent boundary- 
layer flows and, for this reason, the present computer program restricts the 
number of nodes across the boundary layer to 61. Consequently, the chosen 
values of An^ and K must be such that the formula which generates the 
number of grid nodes according to a given or estimated n M , i.e. equation (21), 
does not allow J to exceed 61. Figure 1 is provided, therefore, to provide 
guidance in the selection of J. 



Figure 1. Variation of K with Am for different n -values. 

I oo 


3 



In the present program we march along the t-direction . For a specified 
x-location, the governing equations are solved for each specified t-station. 
Since a linearized form of the equations is being solved, we iterate at each 
t-station until some prescribed convergence criterion is satisfied. For both 
laminar and turbulent flows we use the wall -shear parameter f" as the basis 

W 

for the convergence criterion. For laminar flows the calculation is termin- 
ated when 

IfifJI < 10’ 4 (22a) 

For turbulent flows the error in f" is expressed as a percentage and the 
calculations terminated when 

26f" 

w i 

..fii ~ pfii < 10 (22b) 

0 w w 

In the present method, the flow is laminar at the leading edge, starting 
either as a flat-plate flow or as a stagnation-point flow, and can become 
turbulent at any specified x-station greater than zero. When the transition 
location is not known experimentally, Michel's transition correlation formula 
given in reference 6 is used. According to this formula, transition is pre- 
dicted by the expression 

R e - 1.174[1 + (22,400/R X )]R°- 46 (23) 

tr ax 

0.1 x 10 6 <R < 40 x 10 6 

A 

Previous comparisons with experimental data have shown that this formula works 
well for steady incompressible flows. 

For a given time-dependent external velocity distribution, according to 
equation (4), and geometric configuration, the present computer program deter- 
mines the profiles of f, f 1 , f", b as well as the axial distributions of the 
boundary-layer parameters 6*, e, c^, R { *, R 0 , R x , H. It also computes the 
reduced frequency u>, the phase angles between wall shear and external veloc- 
ity, the phase angle between displacement thickness and external velocity for 
both laminar and turbulent flows, the in-phase and out-of-phase components 
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of an oscillating turbulent flow. The basic details of the comDUtational pro- 
cedure is described in the next two sections. 
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III. DESCRIPTION OF THE COMPUTER PROGRAM 


3.1 Input 

Essentially tne input tc the computer program consists of five types of 
cards. Card 1, shown below, contains the title of the flow problem under 
consideration. 



Load Sheet for Card 1 

The input is punched as an 80-column alphanumeric field. 


Card 2 requires the following information to be specified. The input is 
in 1013 format. 



Load Sheet for Card 2 


NXT Total number of t-stations to be calculated 

NZT Total number of x-stations to be calculated 

NTR x-station where transition begins. For all laminar flows, or flows 

for which transition is to be calculated input 
NTR > NZT. 

IBDY Specifies whether the flow at x = 0 starts as a flat-plate flow or 
as a stagnation-point flow. 

-1 flat-plate flow 
=2 stagnation-point flow 
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ISD Surface distance flag. If surface distance Is calculated the 

pressure gradient parameter P denoted by P2 Is also calculated 
from the external velocity distribution specified In card 5. 

=1 surface distance calculated 
=2 surface distance Input 

IP2 P 2 flag: This option allows the pressure gradient parameter to be 
either Input or calculated from the given external velocity distri- 
bution. If ISD s 1, IP2 must also be equal to 1. 

*1 P 2 calculated 
=2 P 2 Input 

INTR Transition flag 

=1 transition calculated according to (23) 

=2 transition Input 

KPHA Phase angle flag. This flag provides the calculation of phase angle 
between wall shear and external velocity and phase angle between 
displacement thickness and external velocity for both laminar and 
turbulent flows. 

=0 phase angle is not calculated 
»1 phase angle is calculated 

IPOP Flag for in-phase and out-of-phase components of an oscillating 
turbulent flow on a flat plate. 

=0 not computed 
=1 computed 

N Number of t-stations in one cycle to perform the phase angles or 

in-phase angles or out-of-phase components. 

Card 3 requires the following information to be specified. The input is 

in 7F10.0 format. 
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Load Sheet for Card 3 





Load Sheet for Card 3 (continued) 


ETAE Transformed boundary-layer thickness, n B 

DETA(l) Initial An spacing, at^ 

VGP Variable grid parameter, K 

UINF Reference freestream velocity, (ft/sec) 

CB Amplitude of the fluctuating external velocity B 

OMEGA Radial frequency <*> (rad/sec) 

CA Slope of velocity at the stagnation point of an airfoil, (du o /dx) x=0 

Card 4 contains the input t-stations which are read in F10.0 format. 



Load Sheet for Card 4 

Card 5 contains the geometry of the two-dimensional body and the steady- 
state external velocity distribution u Q (x), or the dimensionless pressure 
gradient parameter P(x) in 3F10.0 format. The geometry of the body is 
either read in by specifying the surface distance of the body or is computed 
from the (x/c) and (y/c) coordinates of the body. Here c is a 
reference length of the body, which for an airfoil Is usually taken as the 
chord. The computer names for the usual (x/c), (y/c) coordinates are 
denoted by ZC and YC, respectively. Note that when surface distance 
denoted by Z in field 1, Is to be calculated, then the third field must 
contain the external velocity distribution from which the pressure gradient 
P2 must be calculated. If surface distance Z is Input, and P2 is to be 
calculated, the third field may be blank. There will be NZT cards of this tvoe. 
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L ad Sheet for Card 5 


3.2 Output 

Station header 

NX number of t-s+ation 

NZ number of x-station 

X t, time coordinate 

Z x, space coordinate 

Iteration print 

V(WALL) f; 

DELVW Af^ 

Transition values - If INTR=1 (Card 2), these values will be printed 
close to transition 

RTHETA R 0 , determined from the boundary-layer calculations. 

RTHT (Re) tr , com P ut ed according to formula (23) for a given R v . 

If transition occurs and a station is interpolated, we print out "ew values 
of P2, UO and Z which correspond to P, U Q and x. 

The output boundary-layer parameters include profides f, f', f" and b 
as a function of the similarity variable rj and grid parameter j. Here 

ETA n 

F f ’ ORIGINAL PAGE IS 

U f ' OP POOR QUALITY 

V f" 

B b (=1 + e*) equals 1.0 for laminar flows 

Thev also include disnlacement thickness fi*, momentum thickness e, local 
skin-friction coefficient Cp Reynolds numbers based on 6*, e and x, that 
is, R { *, R q and R^ shape factor H. The definition of these parameters and 
their computer notation is 
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In addition to the above data, the computer program also prints out the 
reduced frequency u (h wx/u 0 ), the phase angle between wall shear and external 
velocity, phase angle between displacement thickness and external velocity for 
both l.iminar and turbulent flows. For turbulent flows it also computes the 
in-phase and ou.-of-phase components of an oscillating turbulent flow on a 
flat platp c . described in section 4.4. 
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IV. SAMPLE CALCULATIONS 


The present computer program can be used to determine the properties of 
two-dimensional laminar and turbulent boundary layers for steady and unsteady 
flows. In order to demonstrate its capability and its input requirements, a 
number of sample calculations are presented in the following subsections. 

4.1 Steady Laminar Flows 

The nonsimilar laminar flow for which the inviscid velocity distribution 
is given by 

u o = u .(' - §■ x) (24) 

is commonly referred to as Howarth's flow and has been computed extensively 
and accurately. Since this is a steady flow, we choose NXT=1 . According to 
previous calculations, flow separation takes place at x = 0.96 and thus, by 
choosing ax to be uniform and equal to 0.05, there are 21 x-stations to 
x = 1.0. Consequently, we set NZT=21 . By setting NTR greater than NZT, say 
50, we avoid the turbulent flow calculations. Since this flow starts as a 
flat-plate flow (P2 = 0), we let IBDY=1. We also set ISD=2, INTR=2 and by 
choosing IRZ=1, we ask the program to compute P2. We also set KPHA, IP0P »N 
equal to zero. 

The transformed boundary- layer thickness is constant for most steady 
laminar flows and a value of 8 is sufficient. For this reason a value of 8 
is read in for ETAE at the first x-station. If needed, for other x-stations, 
the boundary-layer thickness would grow internally. 

In general 40 to 50 node points across the boundary layer are sufficient 
and, for laminar flows, VGP should be set equal to 1.0. Therefore, by choosing 
VGP=1.0 and DETA ( 1 )=0.20, we take 41 points uniformly spaced across the 
boundary layer. We let UINF=1.0. Because this is a steady laminar flow, we 
set CB=0, 0MEGA=0. Also we set CA=0 since the flow starts as a flat-Dlate 
flow. 
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In card 4, we set t( a x)=0, and In card 5, since the surface distance is 
read in, we fill in the 21 x-values as 0, 0.05, 0.10, etc. under z, and the 
values of U0 as 0, 0.9937, 0.98750, etc. under Up. 

Figure 2 shows the computed wall shear parameter f"(0) as a function of 
x. The results agree well with those reported in the literature, see for 
example [8]. The wall shear parameter f“(0), which is positive at x * 0.95, 
becomes negative at the next input station, x - 1.0, indicating that flow sepa- 
ration has taken place between 0.95 < x < 1.0. According to the extrapolation 
shown in figure 1, the flow separation takes place at x - 0.96, and this agrees 
well with other calculations. When f"(0) becomes negative, the calculations 
are terminated internally. 


4.2 Steady Turbulent Flows 

To test the computer program for steady turbulent flows, we have chosen 
three experimental incompressible flows from the data reported at the Stanford 
Conference on Computation of Turbulent Boundary Layers [11], i.e. the flows 
identified as 1400, 2100 and 3300. 



Fiqure 2. Computed wall shear parameter variation with x for Howarth’s flow. 
The circles denote those computed by Cebeci and Smith [8]. 
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Flow 1400 corresponds to a boundary layer with zero pressure gradient and 
the experimental data is due to Wieghart. Flow 2100 corresponds to a boundary 
layer with an initially favorable pressure gradient followed by a zero pressure 
gradient and an adverse pressure gradient that causes the flow to separate; the 
experimental data is due to Schubauer and Klebanoff. Flow 3300 is an equilibrium 
boundary layer with an adverse pressure gradient, and is due to Bradshaw. 

In making these calculations we have first computed a zero pressure gradient 
flow to match the initial conditions provided by the experiment. Then we 
specified the measured external velocity distribution as a function of x and 
performed the boundary- layer calculations. In all cases, we have taken an^ = 0.01 
and K = 1.3. 

For flow 1400, an arbitrary x-distribution was selected with x * 0, 0.01, 
0.05, 0.10, 0.25, 0.50, 0.75, 1, 1.6, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 15, 18, 

21 and u g = 108 ft/sec. The transition location was specified at x = 0.01 
and the calculations performed for NXT=1 and NZT=22. 

Figure 3 allows a comparison of calculated and experimental values of the 
local skin-friction coefficient, c^, as a function of R 0 . The agreement 
between calculated results and experimental results is similar to those shown 
by Cebeci and Smith calculations for steady flows [8]. 

For flow 2100, we have chosen a flat-plate of length approximately 2 ft 
with u g * 117 ft/sec to match the experimental velocity profiles at x « 2 ft. 
Further downstream, the experimental velocity distribution, u e> was 
specified as a function of x and the calculations performed for a total of 
30 x-stations (NZT=30) and for one-time station t = 0 (NXT=1). Figure 4 
allows a comparison of the calculated and experimental shape factor H, Reynolds 
number based on momentum thickness, R„, and local skin-friction coefficient, 

o 

t*. As before, the agreement between calculated and experimental results is 
similar to that demonstrated previously by Cebeci and Smith in steady flow 
calculations [8]. 
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For flow 3300 we have again chosen a flat-plate of length 6 ft and 
tr * 125.63 ft/sec to match the first experimental velocity profile given by 
Bradshaw at x = 2 ft. After that, the experimental velocity distribution was 
specified and the calculations performed for a total of 22 x-stations (NZT=22), 
14 of them corresponding to the flat-plate flow. Figure 5 shows the calculated 
and experimental shape factor H, Reynolds number based on momentum thickness, 
R 0 , and local skin-friction coefficient, c^. 


4.3 Unsteady Laminar Flow on a Flat Plate 


According to Lighthill's analysis [12], the phase anole 4) between the wall 
shear and the external velocity distribution for a laminar flow on a flat plate 
for which the external flow is oscillating according to equation (4) is given 
by separate formulas depending on whether the reduced frequency w(=wx/u 0 ) is 
much smaller or much greater than unity. In order to test the present method 
for this flow and with a range of reduced frequencies, the following procedure 
was adopted to determine the phase angle between u e (x Q ,t) and f^(x Q ,t). The 
values of u Q ( x Q ) and f£"(x ) are first computed from the expressions 


V p 


u_(x ) ■ - f u (x ,t)dt 
oo p J e o 


(31a) 


1 V P 

T X T = p / f ;<v t)dt 


(31b) 


where p * 2tr/w and represents the period of the oscillation of the mainstream. 


From equation (4), 


u e (x Q .t) ~ u o (x o ) = A COSwt (32) 

where A = u 0 (x Q )B. Similarly, we can write 

f;( V t) - f“(x ) = C cos[u)t + <j>(x Q )] = C[coswt cos*(x 0 ) - sinwt si n«t> (x Q ) ] 

( 33 ) 
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with <d(xJ denoting the phase angle between u a and f" at x * x . 

0 C Vi Q 

Integration of the product of equations (32) and (33) leads to 


Vp 


7 


u 0 (x 0 )] • £f“(x 0 .t) 


f;(x 0 )] dt 


cos 4 -(x 0 ) 


ACir/u 


( 34 ) 


with 


V p 


A = “ f [u a (x n t) 
u j u e o» 


U 0 (x„)] dt 


(35a) 


V p 


c -7 7 Cf"(x„,t) -r“Cx 0 )3 z <it 


(35b) 


In order to perform the calculations for one cycle, we must choose At 
according to the formula 

At = ^ _ -j j (36a) 

where N corresponds to the total number of t-statlons, that Is, NXT. Intro- 
ducing the definition of p into equation (36a), we get 

at ■ wr^ry < 36b > 

and compute 

t n * nAt n n ** 1, 2, ..., N (36c) 

In our calculations for this case, we have chosen B * 0.125, u ■ 1.57 rad/sec, 
u Q = 160 ft/sec, An^ = 0.25, K = 1.0 and computed At from (36b). For 
N * 20, corresponding to one cycle, equation (36b) gives 
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At = 1.57(20) = °' 20 


Then t Q =0, t 1 * 0.20, t 2 = 0.40, .... t N = (21 - 1 )At * 4.0. If we want 
to perform the calculations for two cycles, then NXT=41 with At * 0.20 so 
that t^ = 8.0. For three cycles, NXT=61 with At * 0.2 and t^ = 12.0. 

Figure 6 presents results obtained from the present procedure and those 
predicted by Lighthill's analysis: as can be seen, the results are in close 

agreement at the two asymptotes. The res-lts are also in good agreement with 
the numerical computations of Ackerberg and Phillips [13]. They show that if 
B £0.125, Lighthill's low-frequency approximation to the phase angle is satis- 
factory for w < 0.2 and his high-frequency approximation is satisfactory 
for w £ 2.6. 

The calculations _!iown in Table 1 also indicate that it is desir- 
able to compute or two periods in order to determine the phase angle 
accurately. The results of figure 6 were obtained by computing two periods. 
Additional calc. nations indicated that the phase angles computed for three 
periods did not differ from those computed for two periods. 



Fiqure 6. ComDUted phase anqle between wall shear and external velocity for a 
laminar flat-plate flow. 
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Table 1. Computed Phase Angles Between Wall Shear and 
External Velocity According to (4) 


% 

0) 

One 

Cycle 

Two 

Cycles 

0.157 

14.64 

14.64 

0.314 

25.40 

25.30 

0.477 

32.25 

31.79 

0.628 

36.28 

35.20 

0.785 

38.86 

36.99 

C.942 

40.60 

37.91 

1.099 

42.00 

38.58 

1.256 

43.13 

39.20 

1.413 

44.13 

39.95 

1.57 

44.96 

40.77 

1 . J 

45.70 

41.64 

1.884 

46.31 

42.45 

2.n41 

46.85 

43.18 

2.198 

47.29 

43.76 

2.355 

47.67 

44.23 

2.512 

47.98 

44.57 

2.669 

48.26 

44.82 

2.826 

48.49 

44.98 

2.983 

48.70 

45.09 

3.140 

48.87 

45.13 


The procedure used for calculating the phase angle between the wall shear 
and the external velocity distribution can also be used to compute the phase 
angle between the displacement thickness, £*, 

00 

«* =/(' -S-M» 

o e 

and the external velocity. This is achieved by replacing equation (31b) by 

1 V p 

s’TxJJT « ± f 6*(x 0 ,t)dt 

and modifying equations (33), (34) and (35b) in a similar way. In general, at 
low frequencies and at high frequencies, the computed phase angles for one 
cycle and for two cycles are nearly the same (see Table 2). At moderate fre- 
quencies the ones computed by using tw ycles differ from those 
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Table 2. Computed Phase Angles Between Displacement 
Thickness and External Velocity According 
to (4) 


♦j* 


tl> 

One 

Cycle 

Two 

Cycl es 

0.157 

169.56 

169.57 

0.314 

163.03 

163.16 

0.477 

160.06 

160.71 

0.628 

159.70 

161.20 

0.785 

160.73 

163.42 

0.942 

162.35 

166.50 

1.099 

164.12 

1 b'a . 85 

1.256 

165.87 

173.06 

1.413 

167.50 

175.80 

1.57 

168.99 

177.58 

1.727 

170.35 

177.68 

1.884 

171.57 

176.88 

2.041 

172.66 

176.18 

2.198 

173.60 

175.76 

2.355 

174.37 

175.60 

2.512 

174.92 

175.56 

2.669 

175.30 

175.53 

2.826 

175.49 

175.42 

2.983 

175.49 

175.05 

3.140 

175.40 

174.47 


obtained by using one cycle. The difference is probably due to the truncation 
error and can be avoided bv taking more points across the boundary layer. 


Figure 7 presents the computed phase angle results, using one cycle for 
this case. It is seen that at about the same reduced frequency w = 3.0 for 
which reaches it maximum value of 45°, reaches its maximum value of 

175.5° (one cycle). 


4.4 Unsteady Turbulent Flow on a Flat Plate 

An unsteady turbulent boundary-layer flow on a flat plate is considered 
in this subsection which indicates how the computer program can be used to cal- 
culate the flow characteristics corresponding to the configuration of Karlsson 
[14]. The procedure for computing the in-phase and out-of-phase components of 
an oscillating turbulent flow over a flat plate is also outlined. For simplicity. 
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Figure 7. Computed phase angle between displacement thickness and external 
velocity for a laminar flat-plate flow. 

Karlsson's notation is used and denotes the x-component of the velocicy within 
the boundary layer by 

u(x,y,t) * u(x,y) + cos* coswt - u^ sin* sinwt (37) 

Since f' = u/u Q , this expression can be written as 

u Q f' = u(x,y) + u^ cos* coswt - sin<*> sinwt (38) 

Using equation (4) and denoting u B = u^, equation (38) can be written as 


u _ f' 
u^T " B 


°o B 


cos* 


coswt - - — si nut 
u' ’ 


(39) 


To compute the in-phase component cos */u^ from equation (39) we 
multiply both sides of that expression by coswt and integrate with respect to 
t from 0 to to obtain 


u 



co 



f‘(n,t) coswtdt 


(40) 


Similarly, we multiply both sides by sinwt and integrate to determine 
the out-of-phase component, i.e. 
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(41) 



n,t) sinwtdt 


As in the phase-angle calculations, the in-phase and out-of-phase component* of 
the oscillating flow must be computed with uniform steps in At as discussed in 
section 4.3. We approximate equations (40) and (41) by [recalling (f ' )^ 


u™ cos* _ 2 

■ U T1 ) 

oo 

u (1 > sin* _ 2 

-.ii r 


N-l 


X < f j )n «-*» 

(42) 

n=0 


N-l 


T 7 2 (f j )n sinwt n 

(43) 

n=0 



Figure 8 shows the calculated values of the in-phase and out-of-phase components 
for upVu Q = 14. 1 % and u * 25.133 radians/sec or u/2n * 4 cycles/sec. 

These calculations were started as laminar at x = 0 for the external flow 
given by equation (4) with u Q - 17.5 ft/sec. The Ax-spacing for turbulent 
flows was one foot. The turbulent flow calculations were started at x = 0.01 ft 
and the experimental velocity profile matched at x = 12.5 ft (see fig. 6 of 
[7]). Altogether a total of 14 x-stations and 25 t-stations corresponding to 
two cycles were used. Computed results for one cycle and for two cycles show 
that for all practical purposes the results are the same. 


IN-PHASE 



-to.l OUT-OF-PHASE 


ORIGINAL PAGE B 
OF POOR QUALITY 


Figure 8. Computed in-phase ( ) and out-of -phase ( ) components of an 

oscillating turbulent flow for Karlsson's data. 
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Table 3. Computed Phase Ang es Between Wall Shear and 
External Velocity According to Equation (4) 
for Turbulent Flows, w/2ir = 4 cycles/sec, 

B = 14.7*. 



One 

Two 

'V 

Cycle. 

Cycles 

11.49 

21.83 

20.17 

12.93 

22.44 

21.44 

14.36 

22.92 

21.89 

15.80 

23.36 

22.34 

17.23 

23.61 

22.61 

18.67 

23.87 

23.06 


According to the computed results, we also observe tnat the phase angle 
oetween the displacement thickness and external velocity show some sensitivity 
depending on whether one uses one cycle or two cycles as shown in Table 4. 
However, the difference is relatively snuil and can probably be avoided by 
taking more x-stations and more grid points across the boundary layer. 


Table 4. Computed Phase Angles Between Displacement 
Thickness and External Velocity According 
to Eq. (4) for Turbulent Flows, w/2n = 4 
cycles/sec, B = 14.7%. 


V 


b) 

One 

Cyc^e 

Two 

Cycles 

11. 40 

lfia.o? 

161.80 

12.93 

169.57 

155.11 

14.36 

170.13 

163.51 

15.80 

169.93 

170.04 

17.23 

174.08 

173.96 

18.67 

172.23 

171.25 


4.5 Unsteady Laminar and Turbulent Flow on an Airfoil 


To compute a laminar and turbulent flow on an airfoil in which the external 
velocity distribution varies according to eq. (4), we considered NACA 0012 at 
zero angle of attack. The pressure distribution for steady flow is giv^ngjn 


-ss^ 


O* 
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ref. [15]. A total of 25 x-stations and 15 t-stations were used. In the latter 
case a constant At spacing of 0.2 was used. with u = 1.57 rad/sec, An-, - 0.05, 

K * 1.2, u b = 160 ft/sec, and B = 0.125. The slope of u Q (du Q /dx) at the 
stagnation point was determined to be 200/sec. The airfoil geometry was read in 
as (x/c), (y/c) and the surface distance calculations for the given u Q (x) 
distribution were performed internally by setting ISD - 1 . The transition point 
was not specified in the input and the computer program was asked to compute it 
by setting INTR=1. 

Figure 9 shows the computed local skin-friction coefficient and figure 10 
the computed displacement thickness values as a function of x for various values 
of t. We note that, according to equation (23), transition occurs between the 
two input x-stations, x = 0.50027 and x = 0.55027; the interpolated x-station 
that corresponds to the transition station is x = 0.5361. The results in 
these figures show slight wiggles for turbulent flows. This is due to the 
neglect of the finite length of the transition region that exists between a 
laminar and turbulent flow. Although this region is small at low Reynolds num- 
ber, it is not at high Reynolds number and the step change to turbulent flow 
represented by the calculations, lead to oscillations. This can be avoided 
by using an intermittency expression such as that described in [8], 
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Figure 9. Comouted local skin- friction distribution for the NACA 0012 airfoil. 



x 

Figure 10. Computed displacement thickness distribution for the NACA 0012 airfoil 



31 


-fflSSffiS 


V. DISCUSSION AND SUMMARY CONCLUSIONS 


In addition to providing a potential user with a guide to the computer program 
embodying the present calculation method, the results provide information of cur- 
rent and future value. It is clear from the steady-state solutions, that accurate 
results can be obtained with an algebraic eddy-viscosity hypothesis for a wide 
range of pressure gradients including severe adverse gradients which lead to 
separation. The execution time, on a CDC 6600, for a typical calculation (e.g. 
that of figure 2) was of the order of 0.09 sec. 

The unsteady laminar results indicate that calculations have to be completed 
through two periods to provide converged and precise results for the phase angles 
between the wall -shear stress and the freestream velocity and between the displace- 
ment thickness and the freestream velocity. As expected, the displacement-thickness 
phase angle is less amplitude dependent than is the shear-stress phase angle and 
has considerably larger value. Two cycles are also necessary to calculate phase 
angles for turbulent flow which are significantly lower than those for laminar 
flow; the calculated values for shear-stress phase angle agree with available 
measurements within experimental accuracy. The execution time to obtain the 
unsteady turbulent results of figures 6 and 7, was approximately 16 sec and 
suggests that the numerical solution of unsteady, three-dimensional equations 
is a practical and not excessively expensive proposition. 

The airfoil results show that the present procedure may readily be applied 
to practical configurations with calculations from a leading to trailing edge. 
Extensions to include the wake and airfoil separation require the incorporation 
of known techniques ana, like extensions to three-dimensional flow, offer no 
obstacles of principle. The precision of the calculation of transition from 
laminar to turbulent flow is, however, a major unknown. The review of Loehrke, 
Morkovin and Fejer [16] suggests, as might be expected, that the transition 
Reynolds number decreases with the amplitude of freestream fluctuations and is 
also influenced by their frequency. The results upon which conclusions miqht be 
based are not in quantitative agreement but it is orobable that Michel's steady- 
state transition [8] overestimates the transition Reynolds number of up to a 

factor of two. 
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The following paragraphs provide a summary of the more important conclusions 
which may be extracted from the preceding text: 

1. The calculation method can be used, in the form embodied in the 
presently described computer program, to calculate the properties of 
laminar and turbulent boundary layers with pressure gradient and with 

a freestream velocity which fluctuates in time with arbitrary amplitude. 

2. The times, approximately 0.09 sec and 16 sec of CDC 6600 for 
steady and unsteady calculations are modest and suggest that three- 
dimensional calculations can also be performed with acceptable cost. 

3. The calculated phase angles between wall -shear stress and freestream 
velocity increase with amplitude for laminar flow and are greater than 
the almost constant value of approximately 20 degrees obtained for 
turbulent flow. The corresponding phase angles between displacement 
thickness and freestream velocity are around 170 degrees for both 
laminar and turbulent flow. Two cycles are required to obtain converged 
solutions. 

4. The properties of the boundary layer on an airfoil from the stagnation 
line through transition into the downstream turbulent boundary layer 
have also been determined with a varying-amplitude freestream velocity. 
Experimental information is required to confirm the precision of these 
calculations and particularly to indicate the influence of freestream 
fluctuations or transition. The transition correlation of Michel has 
been used for these calculations but alternatives can readily be 
incorporated. 
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VII. COMPUTER LISTING AND A SAMPLE CALCULATION 


In this section we present a listing of the computer programs and 
to illustrate the use of the computer program we present a sample calculation 
for a two-dimensional laminar flow discussed in Section 4.1. Although the 
calculations in that example were done for a total of 21 NZ-stations, we 
only show the first ten stations to illustrate a typical output of our 
computer program. 
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COMMON /BLC 0/ NXT , N ZT , NX , N Z , NTR , NP , ITMAX , I NTR, IBOY , KPHA , I POP , N , 

1 NPT,CNU,ETAE,VGP,A(61) ,ETA(61) ,DETA(61 ) 

COMMON/BLC1/ CA,CB,OMEGA,OMX,REDFR,X(41 ) , Z ( 30) ,UO(30) ,RZ(30) , 

1 PI (30) ,P2(30) ,P3(41,30) ,UE (41 , 30) 

COMMON/PROF/ DELV(61 ) ,F(61 ,41 ,2) ,U(61 ,41 ,2) , V(61 ,41 ,2) ,B(61 , 41 ,2) 

C 

CALL INTIAL 
WRITE (6 ,9000) 

25 WRITE (T>, 91 00) NX,NZ,X(NX) , Z (NZ) 

OMX = OMEGA+X (NX) 

30 IT =0 

I GROW * 0 

60 IT » IT+1 

IF (IT .LE. ITMAX) GO TO 65 
WRITE (6, 2500) 

STOP 

65 IF (NZ .GE. NTR) CALL EDDY 
IF (NZ .GT. 1) GO TO 70 
CALL BCONX 
GO "’O 95 

70 IF (NX .GT. 1) GO TO 90 
CALL BCONZ 
GO TO 95 
90 CALL COEFG 
95 CALL SOLV3 

IF (V (1 ,NX,2) .GE. 0.0) GO TO 61 
IF (NZ .GE. NTR) GO TO 61 
WRITE (6, 9300) 

STOP 

C CHECK FOR CONVERGENCE 

61 IF (NZ .GE. NTR) GO TO 62 
C - LAMINAR FLOW 

IF(ABS(DELV(1)) .GT. 0.001) GO TO 60 
GO TO 100 

C - TURBULENT FLOW 

62 IF (ABS (DELV ( 1 ) / (V( 1 ,NX, 2) +0 »5*DELV ( 1 ) ) ) .GT. 0.02) GO TO 60 
100 IF (NP .EQ. NPT) GO TO 200 

IF (ABS (V (NP ,NX, 2) ) .LE. 0.001) GO TO 200 

IF ( IGROW .EQ. 1) GO TO 200 

IGROW = 1 

LL =2 

CALL OUTPUT (LL) 

GO TO 30 
200 LL = 1 

CALL OUTPUT (LL) 

GO TO 25 

C 

25C0 FORMAT (1 HO, 1 6 X, 25H ITERATIONS EXCEEDED ITMAX) 

9000 F0RMAT(1H1 ,30H** BOUNDARY LAYER CALCULATIONS//) 

9100 FORMAT ( 1 HO , 4I1NX =,I3,5X,4HNZ *,I3,5X,3HX »,F10. 5,5X, 3HZ =,F10.5) 
9300 FORMAT (1 HO, 3 3H** LAMINAR SEPARATION OCCURRED **) 

END 


0F POOR l 
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SUBROUTINE INTIAL 

COMMON /BLCO/ NXT ,NZT,NX,NZ , NTR , NP , ITMAX , INTR, IBDY , KPHA , inP , N , 

1 NPT,CNU,ETAE,VGP,A(61) , ETA (61), BETA (61) 

COMMON /BLC 1 / CA , CB , OMEGA, OMX, REDFR, X ( 4 1 ) ,Z(30) ,UO(30) ,RZ(30) , 

1 PI (30) ,P2 (30), P3 (4 1, 30) ,UE( 4 1,30) 

COMMON/PROF/ DELV(61 ) ,F (61 ,4 1 , 2) , U ( 6 1 , 4 1 , 2) , V( 6 1 ,4 1 , 2) ,B (6 1 , 41 ,2 ) 
DIMENSION TITLE(20) ,ZC(30) ,YC(30) 


c - 

c 

IBDY 

* 

i 

FLAT PLATE 

c 

IBDY 

= 

2 

AIRFOIL 

c 

ISD 

as 

1 

CALCULATED SURFACE DISTANCE 

c 

ISD 

= 

2 

INPUT SURFACE DISTANCE 

c 

IP2 

= 

1 

P2 CALCULATED 

c 

IP2 

= 

2 

P2 READ IN 

c 

INTR 

= 

1 

CALCULATED TRANSITION 

c 

INTR 

= 

2 

INPUT TRANSITION 

c 

KPHA 

= 

0 

DO NOT CALCULATE PHASE ANGLES 

c 

KPHA 

s 

1 

CALCULATE PHASE ANGLES 

c 

IPOP 

= 

C 

DO NOT CALCULATE PHASE COMPONENTS 

c 

IPOP 

= 

1 

CALCULATE PHASE COMPONENTS 

c 

c - 

N 



NUMBER OF X-STATIONS IN ONE ClCLE 


NPT = 6 1 
CNU = 1.6E-04 
ITMAX = 6 

READ(5,8001) TITLE 
WRITE (6, 901 1) TITLE 

READ (5/ 8000) NXT, NZT, NTR, IBDY, ISD, IP 2, INTR, KPHA, IPOP, N 
READ (5, 8100) ETAE , DETA ( 1 ) ,VGP,UINF,CB, OMEGA, CA 
RE AD (5, 8300) (X (I) , 1=1 ,NXT) 

WRI TE ' 6 , 9 2 00 ) NXT , N ZT , NT R , ETAE , DETA ( 1 ) , VGP , CB , OMEGA , UINF , CA 
GO TO (10,30) , ISD 

10 READ (5 , 8400 ) (ZC (I) , YC ( I) ,UO (I) , 1=1 ,NZT) 

WRITE (6,9400) (I,ZC(I) , YC ( I) , 1=1 ,NZT) 

C CALCULATE Z 

Z(1) = 0.0 

UO(1) = UO(1)*UINF 

IF (NZT .EQ. 1) GO TO 60 

SUM1 = 0.C 

DO 20 1=2, NZT 

UO(I) = UO(I)*UINF 

SUM1 = SUM1 +SQRT ( (ZC(I)-ZC(I-I) )**2+(YC(I) -YC (1-1 ) ) **2) 


20 

Z(I) * SUMl 







GO TO 60 






30 

READ\5, 8400) 

(Z 

(I) 

,UO 

(I) 

,P2 (I) ,1=1, NZT) 

50 

IF ( IP2 .EQ. 

2) 

GO 

TO 

55 


60 

IF (IBDY .EQ. 

1) 

P2 

(D 

= 

0.0 


IF (IBDY .EQ. 

2) 

P2 

(1) 

= 

1.0 

55 

DO 90 1=1, NZT 






IF ( I .EQ. 1) 

GO 

i TO 

' 82 




IF ( IP2 .EQ. 

2) 

GO 

TO 

d 1 



IF ( I .EQ. NZT) 

GO 

TO 

70 



Al = (Z(I)-Z(I-I) ) *(Z(I+1)-Z(I-1)) 

A2 = (Z(I)-Z(I-I) ) *(Z (I+1)-Z(I) ) 
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A3 * (Z(I+1)-Z(I))*(Z(I+1)-Z(I-1)) 

DUDS = -(Z(I+U-Z(I))/A1*UO(I-1) + <Z(I+1)-2.0*Z(I) + 

1 Z(I-1))/A2*UO(I) + (Z(I)-Z(I-1))/A3*UO(I+1) 

GO TO 80 

70 A1 = (Z(I-1)-Z(I-2))*(Z(I)-Z(I-2)) 

A2 = (Z(I-1)-Z(I-2))*(Z(I)-Z(I-1)) 

A3 * (Z (I)-Z(I-1))*(Z(I)-Z(I-2)) 

DUDS = (Z(I)-Z(I-1))/A1*UO(I-2)-(Z(I)-Z(I-2) )/A2*UO(I-1) + 

1 (2.C*Z(I)-Z(I-2)-Z(I-1))/A3*UO(I) 

80 P2(I) * Z ( I) /UO (I) *DUDS 

81 REDFR * OMEGA* Z ( I ) /UO { I ) 

GO TO 84 

82 IF(IBDY .EQ. 1) GO TO 31 
REDFR = OMEGA/CA 

84 R2 { I) » UO (I) *Z (I) /CNU 
DO 85 K=1 , NXT 

OMX * OMEGA+X (K) 

UE (K, I) =UO (I) * (1 . 0+CB*COS (OMX) ) 

IF (K .GT. 1 .AND. I .GT. 1) GO TO 83 
IF (K .EQ. 1 .AND. I .GT. 1) GO TO 73 
IF (K .GT. 1) GO TO 71 
IF (IBDY .EQ. 1) P3(K,I) = 0.0 

IF ( IBCY .EQ. 2) P3 (K, I) = ( 1 .0+CB*COS (OMX) ) **2 
GO TO 85 

71 IF ( IBDY .EQ. 1) P3(K,I) = 0.0 

IF (IBDY .EQ. 2) P3{K,I> * (1 .C+CB*COS(OMX) ) **2-REDFR*CB*SIN(OMX) 
GO TO 85 

73 P3 (K, I) =P2 (I) * (1 .0+CB*COS (OMX) ) **2 
GO TO 85 

83 P3(K,I)=P2(I)*(UE(K,I)/UO(I))**2-CB*REDFR*SIN(OMX) 

85 CONTINUE 
90 CONTINUE 

95 WRITE (6, 90 00) (Z (I) ,UO(I) ,P2(I) ,1=1 ,NZT) 

100 DO 110 1=1 ,NZT 

110 PI (I) = 0.5*(P2(I)-ri.0) 

NX = 1 

MZ = 1 

IF ( (VGP-1 .0) .LE. 0.001) GO TO 105 

NP = ALOG((ETAE/DETA(1))* (VGP-1. 0) +1 .0) /ALOG (VGP) + 1.0 
GO TO 112 

105 NP = ETAE/DETA ( 1 '• + 1.0 

112 IF (NP .LE. NPT) GO TO 115 

WRITE (6, 9300) 

STOP 

115 ETA ( 1 ) = 0.0 

DO 120 J=2 ,NPT 
DETA ( J ) = VGP * DETA ( J- 1 ) 

A (J) = 0.5*DETA(J-1) 

120 ETA ( J) = ETA(J-1)+DETA(J-1) 

ETANPQ= 0. 25*ETA (NP) 

STAU15* 1. 5/ETA (NP) 

DO 1 3 C J=1,NP 
ETAB = ETA(J) /ETA(NP) 

ETAB2 = ETAB* *2 
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F ( J , NX , 2 ) = ETANPQ ( *ETAB2 * ( 3 . 0-0 . 5 *ETAB2 ) 

U (J,NX,2) * 0.5*ETAB* (3.0-ETAB2) 

V (J ,NX, 2) = ETAU 1 5 * ( 1 . 0-ETAB2 ) 

B ( J ,NX,2) = 1.C 
130 CONTINUE 
RETURN 

C 

8000 FORMAT ( 1 013) 

8001 FORMAT ( 2 0A4) 

8100 FORMAT ( 8F10. 0) 

830C FORMAT ( FI 0. 0) 

8400 FORMAT ( 3F1 0» 0) 

9C00 FORMAT (///1HG, 46H** EXTERNAL STEADY STATE VELOCITY DISTRIBUTION/ 

1 1 H ,24H** AND PRESSURE GRADIENT/ 1 HO, 4X, 1HZ, 9X, 2IIUO, 8X, 2HP2/ 

2 (1H , 3F10. 5) ) 

9011 FORMAT ( 1II0,20A4) 

9200 FORMAT (///I HO, 12H** CASE DATA/ 1 HO, 3X, 6HNXT =, 13 , 1 4X , 6HNZT =, 

1 I3,14X,6HNTR =,I3/1H ,3X,6HETAE =,E14. 6, 3X, 6HDETA1*, 

2 E14.6,3X, 6HVGP *,E14.6/1H ,3X,6HCB *,E14.6,3X, 

3 6HOMEGA=,E14.6,3X,6HUINF =,E14.6/1H ,3X,fiHCA »,E14.6) 

9300 FORMAT (1HC,37HNP EXCEEDED NPT — PROGRAM TERMINATED) 

9400 FORMAT(//lHC,22H** INPUT BODY GEOMETRY/ 1 HC , 3H MZ,6X, 3HZ/C, 1 IX, 

1 3HY/C/ (1H ,I3,2E14.6)) 

END 
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SUBROUTINE EDDY 

COMMON /BLC 0/ NXT , NZT , NX , N Z , NTR, NP . ITMAX , INTR , IBDY , KPHA , IPOP , N , 

1 NPT , CNU , ETAE , VGP . A ( 6 1 ) , ETA (61), DETA (61) 

COMMON/BLC1/ CA,CB, OMEGA, OMX EOFR,X(41 ) , Z (30) ,UO (30) , RZ (30) , 

1 PI (30) ,P2 (30) ,Pj 1,30) ,UE(41,30) 

COMMON/PROF/ DELV(61 ) ,F(61 ,41 , 2) ,U(61 , 41 ,2) ,V(61 ,41 ,2) ,B (61 ,41 ,2) 


IFLG = 0 

RZ2 = SQRT(RZ (NZ) ) 

RZ216 * RZ2+0.16 
RZ4 = SQRT(RZ2) 

CN * 1 • C 

CRSQV = CN*RZ4*SQRT(ABS (V ( 1 ,NX, 2) ) ) /26.0 
SUM *0.0 
FI * 0.0 
DO 30 J=2,NP 

F2 * U(J,NX,2)M1.0-U(J,NX,2)) 

SUM * SUM+ (F1+F2) *A(J) 

30 FI * F2 

RT = RZ2*SUM 
IF(RT .LE. 425.) GO TO 35 
IF (RT .GT. 6000.) GO TO 38 
XPI * RT/425.-1 • 0 

PI = .55*(1.0-EXP(-.243*SQRT(XPI)-2.98*XPI)) 

GO TO 40 
35 PI = 0.0 
GO TO 40 
38 PI = .55 

40 EDVO = .0168*0. 55/(1. 0+PI))*RZ2*(U(NP, NX, 2)*ETA(NP)-F(NP, NX, 2)) 
J * 1 

50 IF (IFLG .EQ. 1) GO TO 100 
YOA * CRSQV*ETA(J) 

EL =1.0 

IF ( YOA .LT. 4.0) EL = ( 1 . 0-EXP (-YOA) ) **2 
SDVI = RZ2 1 6*ETA ( J) **2 *V( J ,NX, 2) *EL 
IF (EDVI .LT. EDVO) GO TO 200 
IFLG = 1 
100 EDV = EDVO 
GO TO 300 
200 EDV = EDVI 
300 B(J ,NX,2) = 1.0+EDV 
J = J+1 

IF (J .LE. NP) GO TO 50 

RETURN 

END 


ORIGINAL PAGE IS 
OF POOR QUALITY 
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SUBROUTINE ONX 

COMMON/BLCO/ NXT ,NZT,NX,NZ, NTR , NP , ITMAX , INTR , IBDY , KPHA , IPOP , K , 

1 NPT , CNU , ETAE , VGP , A ( 6 1 ) , ETA (61) ,DETA(61) 

COMMON/BLC1/ CA, CB, OMEGA, OMX, REDFR,X (41) , Z(30) ,UO(30) ,RZ(30) , 

1 Pi (30) ,P2 (30) ,P3 (41 ,30) ,UE(41,30) 

COMMON/PROF/ DELV(61) ,F(61,41,2),U(61,41,2),V(61,41,2),B(61,41,2) 
COMMON/BLCC/ Si (61) ,S2(61) , S3 (6 1 ) , S4 (6 1 ) ,S5(61) ,S6(61) , 

1 Rl (61) ,R2(61) #R3(61) 


CEL = 0.0 

IF (NX .EQ. 1) GO TO 30C 
IF ( IBDY .EQ. 1) GO TO 300 
DELX = X ( NX ) —X ( NX— 1 ) 

CEL = 2.0/ (CA*DELX) 

300 U(NP,NX,2) = 1 .0+CB*COS(OMX) 

P22 = -2. 0*P2 (NZ) 

DO 500 J=2,NP 

UB = G.5*(U(J,NX,2)+U(J-1,NX,2) ) 

VB = 0.5*(V(J,NX,2)+V(J-1,NX,2)) 

FVB = 0.5*(F(J, NX, 2) *V ( J, NX, 2) +F (J-1 , NX, 2 ) *V (J-1 , NX, 2 ) ) 

USB = 0.5* (U(J,NX,2) **2+U(J-1 ,NX,2) **2) 

DERBV = (B (J ,NX, 2) *V( J ,NX, 2) -B (J-1 ,NX,2) *V( J- 1 ,NX, 2) ) /DETA ( J-1 ) 

IF (NX .GT. 1) GO TO 400 
R2B = 0.0 
GO TO 450 

400 CUB « 0.5*(U(J, NX-1, 2)+U(J-1, NX-1, 2)) 

CFVB = 0.5*(F(J , NX-1 ,2 ) *V ( J,NX-1 ,2) +F (J-1 , NX-1 , 2 ) *V (J-1 ,NX-1 , 2) ) 
CUSB = C.5*(U(J, NX-1, 2)**2+U(J-1, NX-1, 2)**2) 

CDERBV= (B( J , NX- 1*2) *V(J,NX-1 , 2) -B ( J-1 ,NX-1 , 2) *V( J-1, NX-1 ,2) )/ 

1 DETA ( J- 1 ) 

CR2B = CDERBV+P1 (NZ)*CFVB-P2 (NZ) *CUSB+P3 (NX-1 ,NZ) 

R2B = -CR2B-CEL*CUB 
C 

450 Pi A = A(J) *P1 (NZ) 

51 (J) = B(J,NX,2)+P1A*F(J,NX,2) 

52 (J) = -B (J- 1 ,NX, 2) +P1A*F(J-1 , NX, 2) 

S3 (J) = PI A*V ( J ,NX, 2) 

S4 ( J) = P1A*V (J-1 ,NX, 2) 

S5 (J) = A(J)*(P22*U(J,NX,2) -CEL) 

S6 (J) = A(J) * (P22*U (J-1 , NX, 2) -CEL) 

C 

Rl (J) = F ( J- 1 ,NX, 2) -F (J, NX, 2) +DETA (J-1 ) *UB 
R3 ( J- 1 ) =U ( J- 1 ,NX, 2) -U ( J , NX, 2) +DETA ( J- 1 ) *VB 
R2 ( J) = DETA(J-1)*(R2B-P3(NX,NZ)-(DERBV+P1 (NZ)*FVB- 
1 P2 (NZ ) *USB-CEL*UB) ) 

500 CONTINUE 

Rl (1) = 0.0 
R2(1) = 0.0 
R3(NP)= 0.0 
RETURN 
END 
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c 


c 

c 


c 

c 


SUBROUTINE BCONZ 

COMMON /B LC 0/ NXT , NZT , NX , NZ , NTR, NP , ITMAX , INTR, IBDY , KPHA , IPOP , N # 

1 MPT,CNU,ETAE, VGP,A(61) ,ETA(61) ,DETA(61) 

COMMON/BLC1/ CA, CB, OMEGA, OKX, REDFR,X ( 41 ) , Z ( 30 ) ,UO ( 30) , RZ ( 30) , 

1 PI (30) ,P2 (30) ,P3 (4 1 , 30) ,UE (41 , 3C) 

COMMON/PROF/ DELV(61) ,F(61 ,41 ,2) ,U (61,41,2) ,V(6 1,41,2) ,B (61,4 1,2) 
COMMON/BLCC/ SI (61 ) ,S2 (61 ) .S3 (61 ) ,S4 (61 ) ,S5 (61 ) ,S6 (61) , 

1 R1 (61),R2(61),R3(61) 


U(NP,NX,2) = 1 ,0+CB*COS(OMX) 

BEL = 0.0 

IF (NZ .GT. 1) BEL = 0. 5* (Z (NZ) +Z (NZ-1 ) ) / (Z (NZ) -Z (NZ-1 ) ) 

PIP = PI (NZ)+BEL 

P2P = P2 (NZ) +BEL 

DO 10C J=2,NP 

DEFINITION OF AVERAGED QUANTITIES 

F3 = C.5*(F(J,NX,2)+F(J-1,NX,2) ) 

UB = 0.5*(U(J,NX,2)+U(J-1,NX,2) ) 

VB = C.5*(V(J,NX,2)+V(J-1,NX,2) ) 

FVB = 0.5*(F(J,NX,2) *V( J,NX, 2) +F (J-1 ,NX,2 ) *V (J-1 , NX, 2) ) 

USB = C.5*(U(J,NX,2)**2+U(J-1,NX,2)**2) 

DERBV = (B ( J , NX, 2) *V ( J , NX, 2) -B ( J-1 ,NX,2)*V(J-1 ,NX, 2) ) /DETA (J-1 ) 

IF (NZ .GT. 1) GO TO 10 

CFB = 0.0 

CUB = C.O 

CVB = C.O 

GO TO 20 

10 CFB = 0.5*(F(J,NX, 1) +F (J-1 ,NX, 1 ) ) 

CUB = C.5*(U(J, NX, 1 ) +U ( J— 1 ,NX, 1 ) ) 

CVB = 0.5*(V(J,NX, 1)+V(J-1,NX, 1)) 

CFVB = 0.5*(F(J,MX,1)*V(J,NX,1)+F(J-1,NX,1)*V(J-1,NX,1) ) 

CUSB = 0.5*(U(J,NX, 1)**2+U(J-1 ,NX,1) **2) 

C DERBV = (B ( J , NX, 1 ) *V ( J ,NX, 1 ) -B (J-1 ,NX,1) *V(J-1 ,NX, 1 ) ) /DETA ( J-1 ) 

COEFFICIENTS OF THE DIFFERENCED MOMENTUM EQUATION 

20 SI (J) = B(J,NX,2) /DETA (J-1 ) + (P1P*F (J ,NX, 2) -BEL*CFB) *0.5 

S2 (J) = -B ( J-1 ,NX,2) /DETA (J-1 )+(P1P*F (J-1 ,NX, 2) — BEL*CFB) *0.5 
S3 ( J) = 0.5*(P1P*V(J,NX,2)+BEL*CVB) 

S4(J) = C.5*(P1P*V(J-1 ,NX,2)+BEL*CVB) 

S5 (J) = -P2P*U(J,MX,2) 

S6 (J) = -P2P*U (J-1 ,NX, 2) 


DEFINITIONS OF RJ 

R1 (J) » F ( J- 1 ,NX, 2 ) -F (J,NX,2) +DETA ( J-1 ) *U3 
R3 (J-1 ) =U ( J- 1 ,NX, 2 ) -U ( J,NX, 2 ) +DETA (J-1 ) *VB 
IF (NZ .EQ. 1) GO TO 30 

CLB = CDER3V+P1 (NZ-1)*CFVB-P2(NZ-1)*CUr.B+P3(NX,NZ-1) 
CRB = -P3 (NX, NZ) +BEL* (CFVB-CUSB) -CLB 
GO TO 40 

30 CRB = -P2 (NZ) 

40 R2(J) * CRB- (DERBV+P1P*FVB-P2P*USB-BEL* (CFB*VB-CVB*FB) ) 
100 CONTINUE 

R1 (1) * C.O 
R2(1) = 0.0 
R3(NP)= 0.C 
RETURN 

END ORIGINAL 

OF POOR 


PAGE IS 

quality 
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SUBROUTINE COEFG 

COMMON/BLCO/ NXT , NZT , NX , N2 , NTR, NP , ITMAX , INTR, IBDY , KPHA, IPOP , N , 

1 NPT , CNU , ET AE , VGP , A ( 6 1 ) , ETA (61), DETA (61) 

COMMON/BLC1/ CA ,CB, OMEGA, OMX, REDFR,X(41 ) , Z (3C) ,UO(30) ,RZ(30) , 

1 PI (30) ,P2(30) ,P3 (41 ,30) ,UE(41,30) 

COMMON/PROF/ DELV (61),F(61,41,2) ,U(61,41,2),V(61,41,2),B(61,41,2) 
COMMON/BLCC/ Si (61 ) ,S2 (61 ) ,S3 (61 ) f S4 (61 ) ,S5(61) ,S6(61) , 

1 R1 (61),R2(61) , R3 ( 6 1 ) 


U (NP ,NX, 2) = 1 . 0+CB*COS (OMX) 

DELX = X(NX)-X(NX-1) 

DELZ = Z (NZ) — Z (NZ-1 ) 

ZB = 0.3*(Z(NZ)+Z(NZ-1)) 

UOB = 0.5* (UO (NZ) +UO (NZ-1 ) ) 

CEL = ZB/DELZ 
BEL = ZB/(DELX*UOB) 

CEL2 = G.5*CEL 

P1B = 0.5*(P1 (NZ)+P1 (NZ-1) ) 

P2B = 0.5* (P2 (NZ) +P2 (NZ-1 ) ) 

P3B = 0.25*(P3(NX,NZ)+P3(NX-1,NZ)+P3(NX,NZ-1)+P3(NX-1,NZ-1)) 

P2B 2 = 2. 0*P2B 

P3B4 = 4 . 0*P3B 
BEL 2 = 2 . 0*BEL 

DO 500 J=2 , NP 

FB = C.5*(F(J,NX,2)+F(J-1 ,NX,2) ) 

FVB = 0.5*(F(J,NX,2) *V (J,NX,2)+F ( J-1 , NX, 2 ) *V (J-1 , NX, 2 ) ) 

FB4 = 0. 5*(F(J, NX-1, 2)+F(J-1, NX-1 ,2) ) 

FVJ2 = F(J, NX, 1)*V(J f WX.1)+F(J, NX-1, 1)*V(J, NX-1, 1)+ 

1 F(J, NX-1, 2)*V(J, NX-1, 2) 

FVJ1 = F(J-1,NX,1)*V(J-1, NX, 1)+F(J-1, NX-1, 1)*V(J-1, NX-1, 1)+ 

1 F (J-1 ,NX-1 ,2) *V(J-1 , NX-1 ,2) 

FBI1 = C.25*(F(J ,NX— 1 , 1 ) +F (J-1 , NX-1 , 1 ) +F (J ,NX, 1 ) +F (J-1 ,NX, 1 ) ) 
FVB234= 0.5* (FVJ2+FVJ1 ) 

UB = 0.5* (U ( J,NX, 2) +U (J-1 ,NX, 2) ) 

USB = 0. 5 * (U ( J ,NX, 2 ) **2+U (J-1 ,NX,2)**2) 

UB2 = 0. 5* (U(J ,NX, 1 ) +U (J-1 ,NX, 1 ) ) 

UB4 = 0. 5*(U(J, NX-1, 2)+U(J-1, NX-1, 2) ) 

UJ1 = U(J-1, NX, 1)+U(J-1, NX-1, 1)+U(J-1, NX-1, 2) 

UJ2 = U(J, NX, 1)+U(J, NX-1, 1)+U(J, NX-1, 2) 

UBI1 = 0.25* (U(J ,NX-1 ,1 ) +U (J-1 ,NX— 1 , 1 ) +U (J ,NX, 1 ) +U (J-1 ,NX,1) ) 

UBK1 = 0. 25* (U ( J, NX-1, 1)+U (J-1 ,NX-1 , 1 ) +U (J, NX-1 , 2) +U (J-1 ,NX-1,2) ) 
USJ2 = U(J, NX, 1)**2+U(J, NX-1, 1)**2+U(J, NX-1, 2) **2 
U'JI = U(J-1, NX, 1)**2+U(J-1, NX-1, 1)**2+U(J-1, NX-1, 2)**2 
UB234 = 0.5* (UJ2+UJ1 ) 

USB234* 0 . 5 * ( US J 2 +USJ 1 ) 

VB = 0.5*(V(J, NX, 2 ) +V (J-1 ,NX, 2 ) ) 

VJ1 = V(J-1, NX, 1)+V(J-1, NX-1, 1)+V(J-1, NX-1, 2) 

VJ2 = V(J, NX, 1)+V(J, NX-1, 1)+V(J, NX-1, 2) 

VB234 = 0.5* (VJ2+VJ1 ) 

BVJ1 = B(J-1, NX, 1)*V(J-1, NX, 1)+B(J-1, NX-1, 1)*V(J-1, NX-1, 1) + 

1 B(J-1, NX-1, 2)*V(J-1, NX-1, 2) 

BVJ2 = B(J,NX, 1) *V(J, NX, 1)+B(J, NX-1 ,1 ) *V(J, NX-1 ,1 )+ 

1 B(J, NX-1, 2)*V(J, NX-1, 2) 

DERBV = (B( J ,NX, 2) *V( J ,NX, 2) -B ( J-1 ,NX,2) *V(J-1 ,NX,2) ) /DETA(J-I) 
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CM1 * UB234 

CM2 = UB2-2.0+UBK1 

CM3 * VB234 

CM4 = FB4-2.0*FBI1 

CM5 = FVB234 

CM6 = UB4-2.0*UBI1 

CM7 = CM1*CM6-CM3*CM4 

CM 8 = BVJ2-BVJ1 

CM9 = CM1 +CM6 

CM10 = -CM8+DETA(J-1)*(-P1B*CM5+0SB234*P2B-P3B4+CEL2*CM7+BEL2* 
1 CM2) 

C 

SI (J) = B(J # NX r 2) +A(J) * (P1B*F(J,NX, 2) +CEL2* (FB+CM4) ) 

S2(J) = -B (J-1 ,NX, 2) +A(J) * (PlB*F (J-1 ,NX# 2) +CEL2* (FB+CM4) ) 

S3 ( J) = A(J) * (P1B*V(J ,NX, 2) +CEL2* (VB+CM3) ) 

S4(J) * A(J)*(PlB*V (J-1 ,NX, 2) +CEL2* (VB+CM3) ) 

S5 (J) = A(J ) * (-P2B2*U ( J # NX* 2)— CEL2* ( 2.0*UB+CM9) -BEL2) 

S6 (J) * A(J)M-P2B2*U(J-1,NX,2)-CEL2*(2.0*UB+CM9)-BEL2) 

C 

R1 (J) = F(J-1,NX,2)-F(J,NX,2)+DBTA(J-1)*UB 

R2(J) = CM10-DETA(J-1)*(DERBV+P1B*FVB-P2B*USB-CEL2*(UB*UB+CM9* 

1 UB-VB*FB-CM4*VB-CM3*FB)-BEL2*UB) 

R3 (J-1 )=U (J-1 f NX, 2) -U (J r NX # 2) +DFTA(J-1 ) *VB 
500 CONTINUE 

R3 (NP) = 0.0 
R1 (1 ) = 0.0 
R2 (1 ) = 0.0 
RETURN 
END 


ORIGINAL PAGE IS 
OP POOR QUALITY 
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SUBROUTINE SOLV3 

COMMON/BLCO/ NXT, NZT, NX, NZ , NTR, NP , ITMAX , INTR, IBDY , KPHA , IPOP , N , 

1 NPT , CNU , ETAE , VGP f A ( 6 1 ) ,ETA(61) ,DETA(61 ) 

COMMON/PROF/ DELV (61 ) ,F(61,41,2) ,U(61,41,2) ,V(61,41 , 2 ) , B (61 ,41,2) 
COMMON/BLCC/ SI (61) ,S2(61) ,S3 (61 ) ,S4 (6 1 ) ,S5 (61 ) ,S6 (61) , 

1 R1 (61) ,R2(61) # R3<61) 

DIMENSION All (61 ) , A1 2 (6 1 ) f Al 3 (6 1 ) , A2 1 (61 ) , A22 (61 ) , A23 (61 ) , 

1 G11 (61) ,G12 (61) ,G1 3 ( 61 ) ,G21 (61) ,G22(61) ,G23(61) , 

2 W1 (61) ,W2(61) ,W3(61) ,DELF(61) ,DELU(61) 


W1 (1 ) = R1 (1) 

W2 ( 1 ) = R2 ( 1 ) 

W3 ( 1 ) » R3 ( 1 ) 

All (1)= 1.0 
Al 2 ( 1 ) = 0.0 
Al 3 ( 1 ) = C.O 
A21 (1)= 0.C 
A22(1)= 1.0 
A23 (1 ) - 0.0 
G1 1 (2) =-1.0 
G12 (2)=-0.5*DETA(1 ) 

G1 3 (2) = 0.0 
G21 (2) = S4 (2) 

G23 (2) =-2. 0*S2 (2)/DETA(1 ) 

G22 (2) = G23 (2) +S6 (2) 

DO 500 J=2 , NP 

IF ( J .EQ. 2) GO TO 100 

DEN = (A13(J-1)*A21 ( J-1 ) -A23 (J-1 ) *A1 1 (J-I)-A(J)* 

1 (Al 2 (J-1 ) *A21 (J-1)-A22(J-1)*A11 (J-1) )) 

Gl 1 (J)= (A23 (J-1 ) +A(J) * (A(J) *A2 1 (J-1 ) -A22 ( J-1 ) ) ) /DEN 
G1 2 ( J ) =- ( 1 . 0+G1 1 (J) *Al 1 (J-1) )/A2 1 (J-1) 

Gl 3 (J) = (Gl 1 (J)*A13(J-1)+G12(J) *A23(J-1) )/A(J) 

G21 (J)= (S2(J)*A21 ( J-1 ) -S4 (J) *A23 (J-1)+A(J)*(S4(J)* 

1 A22 (J-D-S6 (J) *A21 (J-1) ))/DEN 

G22 ( J) = (S4 ( J) -G2 1 (J) *Al 1 (J-1) )/A21 (J-1) 

G23 ( J) = (G21 (J)*A1 2 (J-1) +G22 (J) *A22 ( J-1 ) -S6 ( J) ) 

100 All (J) 1.0 

Al 2 (J) =-A(J) -Gl 3 (J) 

Al 3 (J) = A(J) *G1 3 ( J) 

A21 (J)= S3 ( J) 

A22 (J) = S5 ( J) -G23 (J) 

A23 (J) = Si (J) +A(J) *G23 ( J) 

W1 (J) = Rl (J) — Gl 1 ( J) *W1 ( J— 1 ) — Gl 2 ( J) *W2 (J— 1 ) — Gl 3 ( J) *W3 (J— 1 ) 
W2(J) = R2 (J) -G2 1 ( J) *W1 (J-1 ) -G22 (J) *W2 (J-1 ) -G23 (J) *W3 (J-1 ) 
W3 (J) = R3(J) 

500 CONTINUE 


DELU(NP) 

El 

E2 

DELV (NP ) 

1 

DELF (NP) 


W3 (NP) 

W1 (NP) -Al 2 (NP) *DELU (NP) 

W2 (NP) -A22 (NP) *DELU (NP) 

(E2*A1 1 (NP)-EI *A2 1 (NP) ) /(A23 (NP) *A1 1 (NP) -Al 3 (NP) * 
A21 (NP) ) 

(El -Al 3 (NP) * DELV (NP) ) /Al 1 (NP) 
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J * NP 

60C J = J-1 

E3 = W3(J)-DELU(J+1)+A(J+1)*DELV(J+1) 

DELV(J) = (A1 1 ( J) * (W2 ( J) +E3 *A22 ( J) ) -A2 1 ( J) *Wl (J)-E3*A21 (J)*A12(J) 

1 ) / (A2 1 { J) *A1 2 ( J) *A ( J+1 ) -A2 1 ( J) *A1 3 (J) -A(J+1 ) * 

2 A22 (J) *A1 1 (J) +A23 ( J) *A1 1 (J) ) 

DELU(J) =-A(J+1)*DELV(J)-E3 

DELF(J) = (W1 <J)-A12(J)*DELU(J)-A13(J)*DELV(J))/A11 (J) 

IF ( J .GT. 1) GO TO 600 

WRITE (6, 91 00) V(1,WX,2) # DELV(1) 

DO 700 J— 1 /NP 

F(J ,NX,2)~ F (J ,NX, 2) +DELF ( J ) 

U (J ,NX,2) = U ( J , NX, 2 ) +DELU (J) 

700 V (J ,NX,2) = V ( J , NX/ 2 ) +DELV ( J ) 

0(1 /NX/2)= 0.0 
RETURN 



9100 FORMAT (1 II ,5X # 8HV(WALL) =,E14.6,5X,6IIDELVW=/E14.6) 

END 
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SUBROUTINE OUTPUT (LL) 

COMMON/BLCO/ NXT,NZT,NX,NZ ,NTR, WP , ITMAX, INTR, IBDY, KPHA, IPOP ,N , 

1 NPT,CNU,ETAE,VGP,A(61) ,ETA(61) ,DETA(61) 

COMMON/BLC 1 / CA,CB , OMEGA, OMX, REDFR,X ( 4 1 ) , Z ( 30) ,UO(30) ,RZ(30) , 

1 PI (30) ,P2 (30) , P3 (41 ,30) ,UE(41,30) 

COMMON/PROF/ DELV(61) ,F(61,41,2),U(61,41,2),V(61,41,2),B(61,41,2) 
DIMENSION NPK (41) ,RTHT(30) , RTHETA(30) ,ALFAA(61 ,41 ) ,ALFAB(61 ,41 ) , 

1 SUMA (61), SUMB (51) 


IF (LL .EQ. 2) GO TO 150 
NPK (NX) =NP 

IF (RZ (NZ) .EQ. 0.0) GO TO 140 
ZRZ = Z(NZ)/SQRT(RZ(NZ) ) 

DELSTR= ZRZ*(ETA(NP)-F(NP,NX,2)/U(NP,NX,2)) 

CF = 2.0*V(1 ,NX,2)/SQRT(RZ(NZ) ) 

RDELST= UO(NZ) *DELSTR/CNU 
SUM1 = 0.0 

FI = U(1,NX,2)/U(NP,NX,2)*(1.0-U(1,NX,2)/U(NP,NX,2)) 
no 50 J=2 ,NP 

F2 = U(J,NX,2)/U(NP,NX,2)* (1 . 0-U ( J , NX, 2) /U (NP ,NX,2 ) ) 

SUM1 ■ SUM1+ (F1+F2) *A(J) 

50 FI = F2 

THETA = ZRZ*SUM1 

RTIITA = UO(NZ) *THETA/CNU 

H = DELSTR/THETA 

CHECK FOR TRANSITION IF IT IS TO BE CALCULATED 
IF ( INTR .EQ. 2) GO TO 1 50 
IF (P2 (NZ) .GE. 0.0) GO TO 150 
IF (NZ .GE. NTR) GO TO 1 50 
IF (NX .GT. 1 ) GO TO 150 
IF (NZ .EQ. 1) GO TO 150 
IF (NZ .EQ. NZT) GO TO 150 
RZTR = RZ (NZ) 

RTIIETA(NZ) =RTHTA 

RTHT(NZ) = 1 . 174*0 . + 22400. /RZTR) *RZTR**0. 46 
WRITE (6,9 500 ) RTHETA(NZ) , RTHT (NZ) 

IF (RTHETA (NZ) -RTHT (NZ) ) 150,110,120 
110 ZTR = Z (NZ) 

WRITE(6,9600) ZTR 
GO TO 150 

120 ZTR1 = Z (NZ- 1 ) 

ZTR2 = Z (NZ) 

DRTHl = RTHT (NZ-1 ) -RTHETA (NZ-1 ) 

DRTH2 = RTHT (NZ) -RTHETA (NZ) 

ZTR = ZTR1 + (DRTHl * ( ZTR2-ZTR1 ) ) / (DRTHl -DRTH2 ) 

UOTR = UO (NZ- 1 ) + ( (ZTR-ZTR1 ) / (ZTR2-ZTR1 ) ) * (UO (NZ) -UO (NZ-1 ) ) 

P2TR = P2(NZ-1)+( (ZTR-ZTRl ) / (ZTR2-ZTR1 ) ) * (P2 (NZ) -P2 (NZ-1 ) ) 

I = NZT+2 

IF (NZT .EQ. 30) I = NZT+1 
100 I = 1-1 

z (I) = Z ( 1-1 ) 

RZ ( I ) = RZ(I-I) 

UO(I) = UO(I-I) 

PI (I) = PI (1-1) 
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P2 (I) = P2 (I- 1 ) 

DO 90 K«1 ,NXT 
UE(K,I)=UE(K,I-1) 

90 P3(K,I)=P3(K,I-1) 

IF (I .GT. (NZ+1 ) ) GO TO 100 

P2(NZ)= P2TR 

PI (NZ) = 0.5*(1.0+P2(NZ)) 

Z(NZ) = ZTR 

UO (NZ) = UOTR 

RZ(NZ) = Z (NZ) *UO (NZ) /CNU 

UE(1,NZ) = OO(NZ)* (1.0+CB*COS(OMEGA*X(1) ) ) 

P3 (1 ,NZ) = P2(NZ)*(1.0+CB*COS(OMEGA*X(1)))**2 
DO 125 K=2,NXT 

UE(K,NZ) = UO(NZ) *(1.0+CB*COS (OMEGA *X (K) ) ) 

REDFR = OMEGA* Z (NZ) /UO (NZ) 

P3(K,NZ) = P2 (NZ) *(UE(K,NZ)/UO(NZ) ) **2-CB*REDFR*SIN (OMEGA+X (K) ) 
125 CONTINUE 
NTR = NZ 

IF (NZT .LT. 30) NZT = NZT+1 
DO 130 J=1,NP 
F(J,NX,2)= F(J,NX,1) 

U(J,NX,2)= U(J,NX,1) 

V(J,NX,2)= V ( J ,NX, 1 ) 

130 B(J,WX,2)= B(J,NX,1) 

WRITE (6, 9900 ) P2TR,UOTR, ZTR 
IF (NTR .EQ. NZ) RETURN 
140 RTHTA = 0.0 

RTIIETA(NZ) =0. 0 
C 

150 NPO = NP 

NP1 = NP+1 

IF (LL .EQ. 2) NP = NP+2 
IF (NP .GT. NPT) NP * NPT 
DO 160 J=NP1,NPT 

F(J,NX,2) = U (NPO, NX, 2) * (ETA (J) “ETA (NPO) ) +F (NPO, NX, 2) 

U(J,NX,2) = U(NPO,NX,2) 

V (J ,NX,2 ) = V (NPO, NX, 2) 

160 B (J ,NX,2) = B (NPO, NX, 2) 

IF (LL .EQ. 2) RETURN 
C 

180 WRITE (6,9010) 

NPM1 * NP— 1 

WRITE (6, 9000) (J,ETA(J) ,F(J,NX,2) ,U(J,NX,2) ,V(J,NX,2) ,B(J,NX,2) , 
1 J=1 ,NPM1 ,3) 

WRITE (6,9000) NP,ETA (NP) ,F(NP,NX,2) ,U(NP,NX,2) ,V(NP,NX,2) , 

1 B (NP,NX, 2) 

IF (NZ .EQ. 1) GO TO 10 

WRITE (6, 9200) DELSTR, THETA, CF, RDELST, RTHTA, RZ (NZ) ,H, REDFR 
IF (NXT .EQ. 1) GO TO 10 
DELX = X(2)-X(1) 

IF ( IPOP .EQ. 0) GO TO 196 
IF (NZ .LT. NTR) GO TO 1 96 

C CALCULATE IN-PHASE AND OUT-OF-PHASE COMPONENTS OF AN OSCILLATING 
C TURBULENT FLOW 

COMX * COS(OMX) 

SOMX = SIN(OMX) 
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DO 190 J=1,NPT 

ALFAA(J ,NX) = U(J,NX,2)*COMX 

190 ALFAB ( J ,NX) = U ( J ,NX, 2) *SOMX 
IF (NX .LT. NXT) GO TO 10 

11 =2 

12 = N-1 

191 COEFF = 2.0/(CB*FLOAT(N-1)) 

DO 195 J=1 ,NP 

SUMA(J) *0.0 

SUMS (J) =C.C 

DO 198 1*11,12 

SUMA ( J ) =S UMA ( J ) + ALFAA ( J , I ) 

198 SUMS (J) *SUMB (J)+ALFAB( J ,1) 

SUMA(J)= (0.5* (ALFAA(J,I1-1)+ALFAA(J, 12+1) )+SUMA(J) ) *COEFF 

195 SUMS (J) =- (0.5* (ALFAB (J, 11-1 ) +ALFAB ( J, 12+1 ) ) +SUMB (J) ) *COEFF 
WRITE (6, 9300) X(I2+1) , (J,SUMA(J) ,SUMB(J) ,J*1 ,NP) 

IF (12 .EQ. (NXT- 1 ) ) GO TO 196 

11 * 11+ (N— 1) 

12 = 12+ (N-1) 

IF (12 ,LE. (NXT-1 ) ) GO TO 191 
WRITE (6, 9700) 

C 

196 IF (KPHA .EQ. 0) GO TO 10 
IF (NX .LT. NXT) GO TO 10 

C CALCULATE PHASE ANGLES 

11 =2 

12 = N-1 
211 AVI = 0.0 

ADLSTR* 0.0 
DO 210 1=11,12 

ADLSTR= ADLSTR+ (ZRZ* (ETA(NP) -F (NP , 1 , 2 ) /U ( NP , 1 , 2) ) ) 

210 AVI = AVI +V (1,1,2) 

D1 = ZRZ* (ETA(NP) -F(NP,I1-1,2)/U(NP,I1-1,2)) 

D2 = ZRZ* (ETA(NP) -F (NP,I2+1 ,2) /U (NP,I2+1,2) ) 

ADLSTR= OMEGA/6. 2832 *(0. 5 * (D1+D2) +ADLSTR) *DELX 

AVI = OMEGA/6. 2832* (0.5* (V(1 , 11-1 ,2)+V{1, 12+1, 2) ) +AV1 ) *DELX 

ALF2 = 0.0 

BTA2 = 0.0 

ALFBTA= 0.0 

ALFASQ= 0.0 

BETASQ= 0.0 

DO 220 1=11,12 

DLS = ZRZ* (ETA(NP) -F (NP , 1 , 2) /U (NP , 1 , 2 ) ) 

ALF2 = ALF2+ (UE ( I,NZ) -UO (NZ) ) * (DLS-ADLSTR) 

BTA2 = BTA2+ (DLS-ADLSTR)* *2 

ALFBTA= ALFBTA+ (UE ( I,NZ) -UO (NZ) ) * (V ( 1 , I , 2 ) -AVI ) 

ALFASQ* ALFASQ+ (UE(I,NZ) -UO (NZ) ) **2 
220 BETASQ= BETASQ+ (V ( 1 , I , 2) -AVI ) **2 
C CALCULATE PHASE ANGLE BETWEEN WALL SHEAR AND UE 

ALFBTA= (0.5*((UE(I1-1,NZ)-UO(NZ) )* (V ( 1 , I 1 - 1 , 2) -AVI ) + 

1 (UE (12+1 ,NZ ) -UO (NZ) ) *(V( 1,12+1 ,2) -AVI) ) +ALFBTA) *DELX 

ALFASQ* (0.5* ( (UE( 11-1 ,NZ)-UO(NZ) ) **2+(UE(I2+1 ,NZ)-UO(NZ) )**2) + 
1 ALFASQ) *DELX 

BETAS Q= (0.5* ((V(1 , I 1 - 1 ,2 ) -AVI ) ** 2+ (V ( 1 , 12+1 ,2) -AVI) **2)+ 

1 BETASQ) *DELX 

PHI = ARCOS(ALFBTA/SQRT( ALFASQ* BETASQ) ) *57.29578 
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C CALCULATE PHASE ANGLE BETWEEN DISPLACEMENT THICKNESS AND UE 

ALF2 = (0.5*((UE(I1-1,NZ)-UO(NZ))*(Dl-ADLSTR)+(UE(I2+1,NZ)- 
1 UO (NZ) ) * (D2-ADLSTR) ) +ALF2) *DELX 

BTA2 = ( 0.5* ( (D1-ADLSTR) **2+ (D2-ADLSTR) **2) +BTA2) *DELX 
PHI2 = ARCOS (ALF2/SQRT(ALFASQ*BTA2) ) *57.29578 
WRITE (6,9400) X(I2+1 ) ,PHI,PHI2 
C 

IF (12 .EQ. (NXT— 1) ) GO TO 10 

11 = I1+(N-1) 

12 = I2+JN-1) 

IF (12 .LE. (NXT-1) ) GO TO 211 
WRITE (6,9700) 


10 IF (NX .EQ. NXT) GO TO 200 
NX = NX+1 

300 IF (NZ .GT. 1) GO TO 350 
C INITIAL GUESS IN NX DIRECTION (NZ=1) 

310 DO 400 J=1,NPT 

F(J,NX,2)= F(J,NX-1,2) 

U ( J ,NX,2) - U(J,NX-1,2) 

V (J ,NX,2) = V ( J,NX-1 ,2) 

B (J,NX,2) = B ( J , N X- 1 , 2 ) 

400 CONTINUE 
GO TO 370 

350 IF (NX .EQ. 1) GO TO 500 
C INITIAL GUESS IN NX DIRECTION (NZ .GT. 1) 
GO TO 310 

370 IF (NZ .EQ. 1) RETURN 
NP = NPK(NX) 

IF (NX .EQ. 1) RETURN 

IF (NP ,LT. NPK(NX-D) NP = NPK(NX-I) 

RETURN 

200 NX = 1 

WRITE (6, 9800) 

IF (NZ .EQ. NZT) STOP 

NZ = NZ+1 

C SHIFT ALL NX PROFILES IN THE NZ DIRECTION 


DO 

550 

K= 1 

,NXT 



DO 

520 

J=1 

,NPT 



F ( J 

t K 

,1) = 

F ( J i 

,K 

,2) 

U(J 

rK 

,1) = 

U(J, 

,K 

,2) 

V(J 

,K 

rD = 

V(J ( 

. K 

,2) 

B . J 

.K 


B { J i 

rK 

,2) 


520 CONTINUE 
550 CONTINUE 
GO TO 370 

C 

9010 FORMAT (1 HO, 2X, 1HJ , 4X, 3IIETA, 10X,1HF, 13X, 1HU, 13X, 1HV, 13X,1HB) 

9000 FORMAT (1H ,I3,F10.6,4E14. 6) 

9200 FORMAT(1HO, 7HDELSTR=,E14.6, 3X, 7HTHETA *,El 4 . 6 , 3X, 7IICF =, 

1 E14.6/1H , 7HRDELST*,E1 4 . 6, 3X, 7HRTHTA *,E14.6, 3X, 7HRZ =, 

2 E14.6/1H , 7HH =,E14.6,3X, 7HREDFR »,El4.6) 

9300 FORMAT (/1H0 ,4X, 22H** PHASE COMPONENTS **/lH ,18HCYCLE ENDS WITH X* 
1 , E 1 4 . 6/ 1 HO , 2X , 1 HJ , 3X , 8HIN-PHASE , 4X , 1 2HOUT-OF-PHASE/ 

2 ( 1H ,13 ,2E1 4.6)) 
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94C0 FORMAT (1 HO , 1 8HCYCLE ENDS WITH X=,E14.6/ 

1 III , 38HPHASE ANGLE BETWEEN WALL SHEAR AND UE = , 1 2X , El 4 . 6/ 

2 1H , 50IIPHASE ANGLE BETWEEN DISPLACEMENT THICKNESS AND UE= , 

3 E14.6) 

9500 FORMAT (1 HO , 7HRTHETA= , E 1 4 . 6 , 3X , 5HRTHT= , E 1 4 . 6 ) 

9600 FORMAT ( 1 HO ,28H** TRANSITION OCCURRED AT Z=,E14.6) 

9700 FORMAT* 1 HO, 6 7H** INPUT NXT DOES NOT CONTAIN EQUAL INTERVALS OF PHA 
1SE ANGLE PERIOD) 

98C0 FORMAT* 1H0, 39 (2H*-)//) 

9900 FORMAT(1HO,3HP2=,E14.6, 3X , 3HUO=,£l 4 . 6 , 3X,2HZ=,E14. 6) 

END 
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*** TEST CASE 1 - HOWARTH'S LAMINAR FLOW 


** CASE DATA 


NXT = 

1 

N?.T = 

21 

NTR = 

99 

ETAZ = 

0.800000E+01 

DETA1 = 

0.200C0CE+00 

VGP * 

O.ICCOOOE+OI 

CB 

CA 

V 

C.^o0000E+C3 

014EGA* 

C.100000E+01 

UINF = 

O.IGOOOOE+CI 


** EXTERNAL STEADY STATE V7..0CITY DISTRIBUTION 
** AND PRESSURE GRADIENT 


Z 

UO 

P2 

0.0 

1.C0000 

0.0 

0.05000 

0.99375 

-0.00629 

0.10000 

0.98750 

-0.01266 

0.15000 

0.98125 

-0.01911 

0.20000 

0.97500 

-0.02564 

0.2500G 

0.96875 

-0.03226 

0.30000 

0.96250 

-0.03896 

0.35000 

0.95625 

-0.04575 

0.40000 

0.95000 

-0.05263 

0.45000 

0.94375 

-0.05960 

0.50000 

0.93750 

-0.06667 

0.55000 

0.93125 

-0.07369 

0.60000 

0.92500 

-0.08102 

0.65000 

0.91875 

-0.08825 

C. ■'3000 

0.91250 

-0.09579 

0.75000 

0.90625 

-0.10345 

0.80000 

0.90000 

-0.11092 

0.85000 

0.89375 

-0.11880 

C.9C000 

0.88750 

-0.12651 

0.95000 

0.88125 

-0.13461 

1.00000 

0.87500 

-0.14294 

** BOUNDARY 

LAYER CALCULATIONS 
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32 * 


NX 


NX * 

1 NZ 

= 1 X * 

0.0 

Z = 0.0 



V ( WALL ) = 

0.1 87500E+00 

DELVW® G. 

1 98620E+C0 



V ( WALL ) * 

C. 38612CE+00 

DELVW= -0. 

486990E-01 



V (WALL) = 

0. 337421E+00 

DELVW* -C. 

536625E-02 



V (WALL) = 

0. 332055E+C0 

DELVW= -0. 

669 347E-04 


J 

ETA 

F 

U 

V 

B 

1 

0.0 

0.0 

0.0 

G . 331 988E+00 

0. 1C0000E+C1 

4 

0.600000 

0.597017E-C1 

0. 1 98830E+00 

C.3299C3E+C0 

0. 1CC000E+01 

7 

1.200000 

0.237744E+00 

0. 39344 8E+00 

0.31 6348E+0C 

0. ICOOOOE+OI 

10 

1.799999 

0.528925E+C0 

C.574169E+00 

0.282737E+00 

0. 100000E+C1 

13 

2.399999 

0.921 103E+00 

0.728224E+00 

0.228076E+00 

C.100CC0E+01 

16 

2.999998 

0.1 39494E+01 

0.845316E+G0 

0. 161 549E+00 

0.100000E+01 

19 

3.599998 

0.192707E+01 

0.922807E+00 

0.983551E-01 

0. 100000E+01 

22 

4.199997 

0. 24951 9E+01 

0.966684E+00 

0 . 5071 59E-01 

0. 100G00E+C1 

25 

4.799996 

0.308226E+01 

0.987694E+00 

C.219445E-C1 

0.100000E+01 

28 

5.399996 

0.367779E+01 

0.9961 42E+00 

C.792393E-02 

0.100009E+C1 

31 

5.999995 

C.427648E+01 

0.998981E+CC 

0.23791 7E-02 

0. 100C00E+01 

34 

6.599995 

0. 48761 8E+01 

C.999776E+00 

0.592262E-G3 

0. 1000Q0E+01 

37 

7.199994 

0.54761 1E+01 

0. 99996 1E+00 

0.121874E-03 

0.1 0OOOOE+O1 

40 

7.799994 

C.607631E+01 

0.999997E+00 

0 . 206583E-04 

0. 10CCC0E+01 

41 

7.999993 

0.627636E+01 

0.100000E+01 

G.1C9428E-04 

0 . 1 C0C00E+C 1 


ft-*.*-*-*. 


*-*-*-*-*-#-*- 

*-*_*-*-*_*-*- 

*_*_*-*-*-*_* 


NX = 1 NZ = 2 X * 0.0 Z = 0.C50C0 

V (WALL) = C .331 988E+00 DELVW* -0.972310E-C2 

V (WALL) = 0. 322265E+00 DELVW= -0 . 450226E-04 

J ETA F U V B 

1 G.O 0.0 0.0 0 . 322220E+C0 0.100000E+C1 

4 0. 600000 0.581813E-01 0.194094E+00 0.323867E+00 0.100000E+01 

7 1.200000 0.232522E+00 0.386156E+00 0.313812E+00 0.100000E+01 

10 1.799999 0.51 9068E+00 0.566277E+00 0.283169E+00 0.100000E+01 

13 2.399999 0.906733E+00 C.721254E+00 0.23C550E+00 0.100000E+01 

16 2.999998 0 . 1 3 / 691E+C 1 0.840133E+0C 0.164844E+CC 0.100000E+01 

19 3.599998 0.190651E+01 C.919551E+GG C.101347E+00 C . 100000Et-01 

22 4.199997 0.247316E+C1 0.964962E+00 0.527956E-01 C.1C0C00E+01 

25 4.799996 0.305951E+C1 0.986931E+00 G.230901E-C1 0.100000E+G1 

28 5.399996 0.365474E+C1 0.995860E+0Q 0.843160E-02 0.1C0C00E+01 

31 5.999995 0.425332E+01 0.998894E+00 0.256168E-02 0.100000E+01 

34 6.599995 0.485297E+01 0.999754E+00 0.645788E-03 0.100000E+01 

37 7.199994 '.545289E+01 0.999957E+00 0.134743E-03 0.100000E+01 

40 7.799994 0.605288E+01 C.999996E+00 0.232829E-04 0.100000E+01 

41 7.999993 G.625288E+C1 0.100000E+01 0.124532E-04 0.100000E+C1 

DELSTR* 0. 495709E-02 TKETA = 0 . 1 89808E-C2 CF = 0.365695E-01 

RDELST= 0. 307882E+02 RTHTA = 0.117888E+02 RZ = 0.310547E+C3 

= 0.261 164E+01 REDFR = 0.114286E+01 
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NX = 

1 NZ 

= 3 X = 

0.0 

Z = 0.10000 



V (WALL) = 

0.322220E+00 

DELVW* -0. 

998403E-C2 



V (WALL) = 

0.312236E+00 

DELVW= -0. 

657526E-C4 


J 

ETA 

F 

U 

V 

B 

1 

C.O 

0.0 

0.0 

C. 3121 70E+00 

C.100000E+01 

4 

0.600000 

0.5661 33E-01 

0.189204E+00 

0. 31 7602E+00 

C.10000CE+01 

7 

1.200000 

0.227123E+00 

0. 378592E+00 

0.31 1117E+00 

C. 100000E+C1 

10 

1.799999 

0.508849E+00 

0.558047E+00 

0.233535E+00 

0.100000E+01 

13 

2.399999 

0.891792E+00 

C. 71 3937E+00 

C.233062E+C0 

0. 10000CE+C1 

16 

2.999998 

0.135810E+01 

0.834651E+00 

0.168257E+0C 

0. 100000E+01 

19 

3.599998 

0.1 88502E+01 

0.91 6079E+00 

0.104488E+00 

0. 100000E+01 

22 

4.199997 

0.245010E+01 

0.9631 10E+0C 

0.550052E-C1 

0. 100000E+01 

25 

4.799996 

0.303568E+01 

0.9861 03E+00 

0.243212E-C1 

C.10C000E+C1 

28 

5.399996 

0.363058E+01 

0. 995550E+00 

0.898349E-02 

0.1C0000E+01 

31 

5.999995 

0.422905E+01 

0.998798E+C0 

C.276242E-C2 

0.100000E+01 

34 

6.599995 

0.482866E+01 

0.999729E+00 

0.705297E-03 

0. 10000CE+01 

37 

7.199994 

0.542858E+01 

C.999952E+C0 

0.149013E-C3 

0. 100000E+01 

40 

7.799994 

0.602857E+01 

0.999996E+00 

0.261 163E-04 

0. 100000E+01 

41 

7.999993 

0. 622857E+C1 

0. 100000E+01 

0.141 1 65E-C4 

0.100000E+C1 

DELSTR= 0. 71304 0E-02 THETA = 0.271 337E 

-02 CF 

0.251 31 2E-01 

RDELST= 0.440079E+C2 RTIITA = 0.167466E+02 RZ 

0.6171 37E+03 

ll 

= 0.262788E+01 REDFR = 0.114286E+01 



*-*-*-*-*- 




*-*_*-*_*_*_* 


NX = 

1 NZ 

= 4 X = 

o 

• 

o 

Z = 0.15000 



V (WALL) = 

0. 31 21 70E+00 

DELVW= -0. 

102826E-01 



V (WALL) = 

0.301 887E+00 

DELVW= -C. 

81 7980E-04 


J 

ETA 

F 

U 

V 

B 

1 

C.O 

0.0 

0.0 

C.3C1805E+C0 

0. 100000E+01 

4 

0.600000 

0.54991 9E-01 

0.184141E+00 

0.311 074E+C0 

0.1 OOOOOE+OI 

7 

1.200000 

0.221523E+00 

0 . 370720E+00 

0.308233E+00 

0.100000E+01 

10 

1.799999 

0.498217E+00 

0.549423E+00 

0.28331 8E+00 

0.100000E+C1 

13 

2.399999 

0.8761 95E+G0 

0.70621 3E+00 

C.235609E+00 

C. 1 C0000E+01 

16 

2.999998 

0.13384 1J+01 

0.828816E+00 

0.171 804E+00 

0.1 00000E+01 

19 

3.599998 

0.1 86247E+01 

0.91 2348E+00 

0. 107804E+00 

0,1 00000E+C1 

22 

4.199997 

0.242585E+01 

0.961 1 00E+00 

0.573694E-01 

0 . 1 00000E+01 

25 

4.799996 

0. 301 058E+01 

0. 9851 94E+C0 

0.256553E-01 

0. 100000E+01 

28 

5.399996 

0.36051 3E+01 

0.995206E+00 

0.958888E-02 

0.100000E+01 

31 

5.999995 

0.420347E+C1 

0.998690E+00 

0 . 298533E-02 

0.1 00000E+01 

34 

6.599995 

0.480205E+01 

0.999701E+00 

0 . 772255E-03 

0. 100000E+01 

37 

7.199994 

0. 540296E+01 

0.999946E+00 

0.165874E-03 

0.1 00000E+0 1 

40 

7.799994 

0 • 600294E+0 1 

0.999995E+00 

0. 29841 9E-04 

0. 100000E+01 

41 

7.999993 

0.620294E+01 

0. 100000E+01 

C.162542E-C4 

0. 100000E+01 

OELSTR= 0. 88874 3E-02 THETA = 0.335972E 

-02 CF 

0. 19901 3E-01 

RDELST= C.545049E+02 RTHTA = 0.206C45E+02 RZ 

0.91 9922E+03 

H 

= 0.264529E+01 REDFR = 0.114286E+01 







*_*•*•*.♦_*.* 
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;jx = 

1 NZ 

* 5 X = 

0.0 

Z « 0.2C000 



V (WALL) = 

0.301 8O5E+0G 

DELVW= -C. 

1C6050E-01 



V (WALL) = 

0.291200E+00 

DELVW= -C. 

978548E-C4 


j 

ETA 

F 

U 

V 

B 

1 

0.0 

0.0 

0.0 

C.291 103E+00 

0.1000C0E+C1 

4 

0.60C000 

0. 5331 29E-0 1 

0.1 78890E+C0 

0. 304263E+00 

0.100000E+C1 

7 

1.200000 

C.21 5707E+00 

0.362511E+00 

0.3051 43E+C0 

C.100000E+01 

10 

1.799999 

0.4871 34E+00 

0.540372E+00 

0.284004E+00 

0.100000E+C1 

13 

2.399999 

0.859882E+00 

0.69804 1 E+00 

0.2381 90E+00 

0.10000CE+C1 

16 

2.999998 

0.131775E+01 

0.822587E+00 

0. 1 75493E+00 

0. 100000E+C1 

19 

3.599998 

0. 183873E4-01 

C.908327E+C0 

0.11 1 31 0E+G0 

C.10GC00E4-C1 

22 

4.199997 

0.240028E+01 

0.95891 1E+00 

0.599049E-01 

0. 100000E+01 

25 

4.799996 

0.298408E+01 

0.984194E+00 

0.271057E-C1 

C.100000E4-C1 

28 

5.399996 

0.357824E4-01 

0.994824E+00 

0.102561E-01 

0.100000E+01 

31 

5.999995 

0.41 7644E+01 

C.998568E+00 

0.323399E-C2 

0.1 OOOOOE4-C1 

34 

6.599995 

0.477597E4-01 

0.999669E+00 

0.848130E-03 

0.1C0000E+01 

37 

7.199994 

0.537587E+01 

0.999939E+00 

0.1 84847E-03 

0 . 1000C0E+C1 

40 

7.799994 

0. 597586E+01 

0.999995E4-00 

0.337430E-C4 

0.100000E4-01 

41 

7.999993 

0.61 7586E+01 

0. 100000E+01 

0.186796E-C4 

0.1 OOOOOE+OI 


DELST R= 0.104503E-01 THETA = 0.392278E-02 CF = 0.166770E-01 
RDELST= C.636816E+02 RTHTA = 0 .2 39044E4-02 RZ = 0.121875E+C4 
H = 0.266401E+01 REDFR = 0.114286E+01 








.MX = 

1 NZ 

= 6 X = 

0.0 

Z = 0.25000 



V (WALL) = 

0.291 103E+00 

DELVW= -0. 

1 09539E-C1 



V (WALL) = 

0.28014 9E+00 

DELVW= -0. 

1 1 5092E-03 


J 

ETA 

F 

U 

V 

B 

1 

O.C 

0.0 

0.0 

0.28C033E4-00 

0.100000E4-01 

4 

0.600000 

0.51571 3E-01 

0.1 73435E4-00 

0. 297141 E4- 00 

0 * 1 00000E4-0 1 

7 

1.2C0000 

0.209655E4-00 

0.35393524-00 

0. 301 824E4-C0 

0.100000E4-01 

10 

1.799999 

0.475562E4-00 

0.530849E4-00 

0.284080E4-00 

0.100000E4-01 

13 

2.399999 

0.842785E4-C0 

C.689375E4-00 

0.240801E4-00 

0.100000E4-01 

16 

2.999998 

0.129602E4-01 

0.815920E4-00 

0.1 79335E4-00 

0.100000E4-01 

19 

3.599998 

0.181370E4-01 

0.903981E4-00 

0.1 15025E4-00 

0.100000E4-01 

22 

4.199997 

0.237325E4-01 

0*95651 9E4-00 

0.626322E-01 

0.100000E4-01 

25 

4.799996 

0.295604E4-01 

0.983088E4-00 

0.286873E-01 

0.100000E4-01 

28 

5.399996 

0.354976E+01 

0.994396E4-00 

0. 109933 E- Cl 

0.100000E4-01 

31 

5.999995 

0.41 4780E4-01 

0.998430E4-00 

0.35 1 333E-02 

0.100000E4-01 

34 

6.599995 

0.474 72 9E4-01 

0.99963224-00 

0.93421 4E-03 

0. 100000E+01 

37 

7.199994 

C. 53471 8E4-01 

0.999932E4-00 

0 . 206550E-03 

0.100000E4-01 

40 

7.799994 

0.59471 6E4-01 

0.999994E4-00 

0.387246E-04 

0.1 OOOOOE4-01 

41 

7.999993 

0.61 471 6E4-01 

0.100000E4-01 

0.21 4701E-04 

0.100000E4-01 

DELST R= 0.1 19059E-01 THETA = 0.443553E 

-02 CF 

0 . 1 43954E-0 1 

RDELST= 0. 720863E+02 RTHTA = 0.268557E+02 RZ 

0.151367E4-04 

!i 

= 0.268420E+01 REDFR = 0.114286E+01 
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NX * 

1 NZ 

* 7 X = 

0.0 

Z * 0.30000 



V (WALL) * 

0.28C033E+00 

DELVW* -0. 

1 1 3329E-01 



V (WALL) « 

0.2687C1E+00 

DELVW* -0. 

133248E-03 


J 

ETA 

F 

U 

V 

B 

1 

0.0 

0.0 

C.O 

0.268567E+00 

0.100000E+C1 

4 

0.600000 

0.497615E-01 

0.1 67759E+00 

0.289678E+00 

C.100C00E+01 

7 

1.20CC00 

0.203345E+00 

0.. *956E+00 

C.298249E+00 

C.100000E+C1 

10 

1.799999 

0.463452E+00 

0.520807E+O0 

C.284027E+00 

0.100000E+01 

13 

2.399999 

0.824825E+00 

C .6801 58E+00 

0.243439E+00 

0. 1000C0E+01 

16 

2.999998 

0.12731 1E+01 

0.808762E+00 

0.18334CE+00 

0.100000E+01 

19 

3.599998 

0.1 7872 4E+01 

0.899264E+00 

0.11 8970E+00 

0. 100000E+C1 

22 

4.199997 

0.234461E+01 

0.953895E+00 

0.655742E-01 

0.100000E+01 

25 

4.799996 

0. 29262 7E+01 

0.981 861 £9-00 

0.304194E-01 

0.100000E+C1 

28 

5.399996 

0. 351 950E+01 

0.993916E+00 

0.11 81 26E-01 

0.100000E+01 

31 

5.999995 

0.41 1737E+01 

0.998273E+00 

0. 382744E-02 

0. 100000E+01 

34 

6.599995 

0.471680E+01 

0.999590E+00 

0.103287E-02 

0.100000E+C1 

37 

7.199994 

0.531667E+01 

0.999922E+00 

C.232568E-03 

0.100000E+01 

40 

7.799994 

0. 591 665E+01 

0.999993E+00 

G.444706E-04 

0. 100000E+01 

41 

7. 299993 

0.61 1665E+01 

0. 100000E+01 

0.252692E-04 

C. 100000E+01 


DELSTR= 0.1 32999E-01 THETA = 0.491485E-02 CF = 0.126439E-01 

RDELST= 0.800073E+02 RTHTA => 0.295659E+02 RZ = 0.180469E+C4 

H = 0.270607E+01 REDFR = 0.114286E+01 






*-♦-*_*_*.♦_*_ 

* 

1 

* 

1 

* 

1 

* 

t 

* 

1 

# 

1 

* 

NX = 

1 NZ 

= 8 X = 

0.0 

Z = 0.35000 



V (WALL) = 

0.268567E+00 

DELVW= -0. 

1 1 7447E-01 



V (WALL) = 

0.256823E+00 

DELVW= -0. 

155162E-03 


J 

ETA 

F 

U 

V 

B 

1 

0.0 

0.0 

0.0 

0.256667E+00 

0. 100000E+C1 

4 

0.6C0000 

0. 478772E-01 

0.161 839E+00 

0.28184CE+00 

C . 1 OOCOOE+OI 

7 

1.200000 

0. 196752E+00 

0. 335533E+00 

0.294388E+00 

C.100000E+01 

10 

1.799999 

0.450749E+00 

0.5101 89E+00 

0 .283825E+00 

0. 100C00E+01 

13 

2.399999 

0. 805910E+00 

0.670325E+00 

0.246099E+00 

0. 100000E+01 

16 

2.999998 

0.124889E+01 

0.801 049E+00 

0.1 8751 9E+00 

0.100000E+01 

19 

3.599998 

0.175916E+01 

0.894127E+00 

0 . 1 231 72E+00 

0. 100000E+01 

22 

4.199997 

0.231415E+01 

0.951002E+00 

0.687604E-01 

0.100000E+01 

25 

4.799996 

0.289457E+01 

0.980491E+00 

0.323235E-01 

0. 100000E+01 

28 

5.399996 

0.348725E+01 

0.993372E+00 

0.1 27264E-0 1 

0.100000E+01 

31 

5.999995 

0.408492E+01 

0. 99809 3E+00 

0.41 8363E-02 

0. 100000E+01 

34 

6.599995 

0.468429E+01 

0.999540E+00 

0.11461 0E-02 

0. 100000E+01 

37 

7.199994 

0.52841 4E+01 

0.999912E+C0 

0.262093E-03 

0.100000E+01 

40 

7.799994 

0.58841 2E+01 

0.999992E+00 

0.51431 3E-04 

0.100000E+01 

41 

7.999993 

0.60841 2E+01 

0.100000E+01 

0.297837E-G4 

C . 1 000C0E+01 

DELST R= 0.1 4661 4E-01 THETA = 0.537079E 

-02 CF 

0.11 22 38E-01 

RDELST= C.876247E+02 RTHTA = 0.32C989E+02 RZ 

0.2091 80E+C4 

II 

= 0.272984E+01 REDFR = C.114286E+01 







#-.*_* — *—*—#-* 


57 



MX » 

1 NZ 

= 9 X = 

0.0 

Z = C.4C0C0 



V (WALL) = 

0.256667E+00 

DELVW= -0. 

121982E-C1 



V (WALL) = 

0.244469E+00 

DELVW= -0. 

181 634E-03 


J 

ETA 

F 

U 

V 

D 

1 

0.0 

C.O 

0.0 

0.244287E+00 

C.1C0000E+01 

4 

0.600000 

0.459102E-01 

0.1 55649E+00 

0.273584E+00 

0.100000E+01 

7 

1.200000 

0.1 89844E+00 

0.325616E+C0 

C.290204E+C0 

0 • 1 OCOOOE+CI 

10 

1.799999 

0.437384E+00 

0.498923E+0C 

0.283448E+00 

0*1 C0000E+0 1 

13 

2.399999 

0. 785926E+00 

0.659796E+00 

0.248771E+00 

0 . 1 COOOOE+O 1 

16 

2.999998 

C. I22320E+01 

0.792702E+00 

0.191 889E+0C 

C . 1 C0000E L C 1 

19 

3.599998 

C.172927E+01 

0.883503E+00 

0.127658E+C0 

0. 10C00GE+C1 

22 

4.199997 

C.228164E+01 

C.947797E+CC 

0.722225E-C1 

0.100000E+01 

25 

4.799996 

0.286068E+C1 

0.978954E+00 

0.344275E-C1 

C.100000E+G1 

28 

5.399996 

0. 345274E+01 

0.992754E+00 

0.1 37529E-01 

0.1 0000CE+0 1 

31 

5.999995 

0. 40501 7E+01 

0.997885E+0C 

C.458946E-02 

0. 1COOOOE+C1 

34 

6.599995 

0.464947E+01 

0.999483E+0G 

C.127798E-02 

C . 1 C0000E+01 

37 

7.199994 

0. 524931E+C1 

0.999899E+00 

0.297499E-C3 

0. 1 C000CE+G1 

40 

7.799994 

0.584928E+01 

0.99999 1E+00 

0.59681 6E-04 

0. 1G00C0E+01 

41 

7.999993 

0. 604928E+01 

0. 100000E+G1 

0. 350289E-04 

0 . 1 OOOOOE+C 1 

DELSTR= 0.160111 E-01 THETA = 0.580995E 

-02 CF 

0.100253 E-01 

RDKLST= 0. 950 

658E+02 RTHTA = C.344966E+C2 RZ 

0.237500E+C4 

II 

= 0.275581E+01 REDFR = 0.114286E+01 






*-*-*-*-*-*-*- 

*-*-*_*-*-♦-* 


NX = 

1 NZ 

= 10 X = 

0.0 

Z = 0.45000 



V (WALL) = 

0.244287E+00 

DELVW= -0. 

126964E-01 



V (WALL) - 

0.231591E+00 

DELVW= -0. 

212776E-C3 


J 

ETA 

F 

U 

V 

B 

1 

0.0 

0.0 

C.O 

0.231 378E+00 

0.100000E+01 

4 

0.600000 

0.438514E-01 

0.149158E+00 

0.264859E+00 

0.1 00000E+01 

7 

1.200000 

0.1 82585E+00 

0.3151 44E+00 

0.285651 E+00 

0* 1 00000E+C 1 

10 

1.799999 

0. 423279E+00 

0 . 486928E+00 

0.282863E+00 

0.100000E+01 

13 

2.399999 

0. 764739E+00 

0.648475E+00 

0.251448E+00 

0. 1 000C0E+0 1 

16 

2.999998 

0.1 19584E+01 

0.783629E+00 

0. 196463E+00 

0. 1 OOOOOE+OI 

19 

3.599998 

0. 1 69733E+01 

0.88231 5E+00 

C. 1 32466E 4 00 

0. 1 00000£+0 1 

22 

4.199997 

0.224680E+01 

0.944224E+00 

0.760038E-01 

0.1 00000E+0 1 

25 

4.799996 

0.282428E+01 

0 .9772 1 6E+00 

0 . 367653E-01 

0. 100000E+01 

28 

5.399996 

0.341564E+01 

0.992045E+00 

0.149122E-01 

0. 100000E+C1 

31 

5.999995 

0.401 281E+01 

0.997643E+00 

0.505615E-02 

0 . 100000E+0 1 

34 

6.599995 

0.461 202E+0 1 

0.99941 4E+00 

0.1431 37E-02 

0.100000E+'; • 

37 

7.199994 

0. 52 1 1 84E+01 

0 • 999883E+00 

C.338957E-03 

0. 100000E+0 ’ 

40 

7.799994 

0. 581181 E+ 01 

0. 999989E+00 

0.702 64 8E-04 

0.100000E+01 

41 

7.999993 

0.601 181E+01 

0. 100000E+01 

0.411 659E-04 

0. 1000COE+01 

OELST R= 0. 1 73658E-01 THETA = 0.623697E 

-02 CF 

0 .8982 10E-02 

RDELST= 0. 102431 E+ 03 RTHTA = 0.367883E+02 RZ 

0.265430E+04 

H 

= 0. 278434E+01 REDFR = 0.114286E+01 
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1 6 Abstract 


There are many unsteady boundary- layer problems in which it is necessary to account for ] 

fluctuations in the external flow. These fluctuations may change in direction and magnitude and, 
in most turbulent flow cases, may be regarded as low frequency fluctuations superimposed on the 
turbulence energy spectrum. The computer program described here provides solutions of two- 
dimensional equations appropriate to laminar and turbulent boundary layers for boundary conditions 
with an external flow which fluctuates in magnitude; it can readily be extended to the more 
general case. It is based on the numerical solution of the governing boundary-layer equations 
by an efficient two-point finite-difference method developed by Keller and Cebeci [1]. An eddy- 
viscosity formulation is used to model the Reynolds shear-stress term. Here, after briefly 
describing the main feature r of the method, we provide instructions for the computer orogram with 
a listing and present sample calculations to demonstrate its usage and capabilities for laminar 
and turbulent unsteady boundary layers with an extern 1 flow which fluctuates in magnitude. For 
further details of the method, the reader is referred i reference 2. \ 


-SMS# 


17 Key Words (Suggested by Author(s)) 

Unsteady boundary layer 

Laminar and turbulent boundary layer 

Boundary-layer computation 

18. Distribution Statement 

Uni imi ted 

STAR Category - 34 


19 Security Oassif (of this report) 

20 Security Oassif. (of this page) 

21 No of Pages 

22 Price* 

Unclassified 

Unclassified 


61 

S4.25 


'For sale by the National Technical Information Service. Springfield. Virgtnn 22161 





















